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ABSTRACT 


ALGEBRAIC GRID GENERATION 
USING 

TENSOR PRODUCT B-SPLINES 

Bonita Valerie Saunders 
Old Dominion University, 1985 
Director: Dr. Philip W. Smith 

In general, finite difference methods are more success- 
ful if the accompanying grid has lines which are smooth 
and nearly orthogonal. This thesis discusses the develop- 
ment of an algorithm which produces such a grid when given 
the boundary description. 

Topological considerations in structuring the grid 
generation mapping are discussed. In particular, this 
thesis examines the concept of the degree of a mapping 
and how it can be used to determine what requirements are 
necessary if a mapping is to produce a suitable grid. 

The grid generation algorithm uses a mapping composed 
of bicubic B-splines. Bounaary coefficients are chosen 
so that the splines produce Schoenberg's variation diminish- 
ing spline approximation to the boundary. Interior coeffi- 
cients are initially chosen to give a variation diminishing 
approx i mat i on to the transfinite bilinear interpolant of 
the function mapping the boundary of the unit square onto 


i 


the boundary of the grid. 

The practicality of optimizing the grid by minimizing 
a functional involving the Jacobian of the grid generation 
mapping at each interior grid point and the dot product 
of vectors tangent to the grid lines is investigated. 

Grids generated by using the algorithm are presented. 
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1. INTRODUCTION 


Grid generation is the numerical development of curvi- 
linear coordinate systems. In recent years grid generation 
has been the key to solving partial differential equations 
on arbitrarily shaped regions by finite difference methods. 
Although much of the motivation for grid generation has 
come from fluid dynamics, the techniques apply to any area, 
such as electromagnetics and heat transfer, which involves 
the solving of partial differential equations on a physical 
domain. 

Inherent in grid generation techniques is a mapping 
T from some canonical domain such as a square or rectangle 
in two dimensions, or cube in three dimensions, onto the 
physical domain on which the partial differential equations 
are to be solved. The image of a mesh on the canonical, 
or computational, domain will be a grid on the physical 
domain. When the grid boundary coincides with the boundary 
of the physical domain, the system generated is called 
a boundary fitted coordinate system. 

A boundary fitted coordinate system allows one to 
apply boundary conditions exactly, thus avoiding interpola- 
tion errors. However, such a system may make the equations 
to be solved more complex [Sm]. 

The distribution of the coordinate lines, or grid 
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lines, should be smooth, but concentrated in areas where a 
large gradient occurs in the physical solution. As stated 
by Thompson, Warsi and Mastin [ TWM] , "the grid points may 
be thought of as a finite set of observers of the physical 
solution, stationed to be most effective in covering all 
of the action on the field." Ideally, the grid should 
be adaptive, that is, coupled with the physical solution 
so that it automatically redistributes its grid lines to 
obtain the desired regions of concentration as the solution 
evolves. However, the interior lines should not cross 
the physical boundary and should be nearly orthogonal at 
the intersection points to avoid large truncation errors 
in the finite difference approximations. 

Grid generation is based on the observation that 
finite difference computations are much easier to make 
on a uniform mesh over a canonical domain such as a square 
or cube than on a grid over an irregularly shaped region. 
Therefore, the partial differential equations to be solved 
must first be transformed so that the computational coordi- 
nates become the independent coordinates. The resulting 
equations may then be Expressed as finite difference equa- 
tions on the computational domain. 

Grid generation techniques may be divided into two 
general types: partial differential equation methods and 

algebraic methods. P.d.e. methods include elliptic, hyper- 
bolic and conformal mapping techniques. All of these methods 
involve the solving of partial differential equations to 


obtain the grid coordinates. The simplest elliptic method 
for grid generation uses the Laplace equations 
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A* C 


a* n 



0 

0 


where c and n are The computational coordinates and x and y 
are the physical coordinates in two dimensions. The equations 
are first transformed so that the independent and dependent 
variables are interchanged. Then the new equations are 
solved for x and y in terms of e and n. Some control over 
the grid cell spacing can be accomplished by introducing 
control functions P(e,n). QU.n) and solving the Poisson 
equations [TWM, p. 39] 

A a e » PU.n) 

A* n 3 Q ( e . n ) . 

Solving the Laplace equations 

a* C * 0 

A* n » 0 

with boundary conditions 

» 

C X 3 Hy 

h * -n x 

produces a conformal transformation [TWM, p. 11] 

Starius [St, p. 27] shows that solving an initial value 
problem satisfying 
x n s -n F 
yn « x e F 
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where F is chosen so that the system is hyperbolic produces 
a hyperbolic grid generating system. Grids generated from 
elliptic equations are generally smooth regardless of the 
type of boundary, but slope discontinuities propagate through 
hyperbol ical ly generated grids [5tj. Generating a grid 
using conformal mapping techniques requires careful selec- 
tion of the boundary data, making it difficult to structure 
the grid to obtain a high concentration of grid points 
in areas of large gradients in the physical solution. More 
grid points may have to be added in order to capture regions 
of rapid change such as shocks and boundary layers. Also 
in p.d.e. generated systems the Jacobian information needed 
for the transf ormat i on of the equations being solved must 
be computed numerically. 

in algebraic methods an explicit functional relation- 
ship between the computational and physical domains is 
defined. Therefore, no p.d.e. need be solved to obtain 
the grid coordinates and the Jacobian matrix can be computed 
analytically. Such methods allow more precise controls 
of the grid structure making it easier to concentrate grid 
points in large gradient areas. However, a 1 gebra 1 ca 1 1 y 
generated grids are more sensitive to point distributions 
on the boundary and, in general, may not be as smooth as 
those generated by elliptic techniques [Sm]. Slope discon- 
tinuities on the boundary may propagate into the field. 
Nevertheless, a variety of techniques have been used to 
produce acceptable smoothness in algebraically generated grids. 
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TMs thesis discusses an algebraic grid generation 
technique for creating boundary fitted coordinate systems. 
This technique uses a mapping which is a sum of tensor 
product B-splines. Chapter 2 discusses degree theory, 
explaining hew the degree of a mapping can be used to deter- 
mine what conditions must be met if an algebraic transforma- 
tion is to produce a suitable grid. Chapter 3 presents 
the tensor product grid generation mapping and discusses 
the properties of B-splines to show their suitability for 
use in such a mapping. Chapter 3 also introduces a func- 
tional which can be used to change the coefficients in 
the mapping in order to enhance the smoothness and orthogo- 
nality in the generated grid. 

Chapter 4 discusses the computer program TENTEST 
which uses the techniques p-^sented in Chapter 3 to generate 
grids on erbitrarily shaped two-dimensional domains. Some 
of the grids created using TENTEST are illustrated and 
discussed in Chapter 5. Conclusions and suggestions for 
further study are presented in Chapter 6. 


2. APPLICATIONS OF DEGREE THEORY 


This chapter discusses degree thoery and shows how 
the degree of a mapping can be used to help determine what 
requirements are necessary if a transformation T is to 
produce a suitable grid. 

Since the distribution of grid lines should be smooth 
with concentration in areas of large gradients in the 
physical solution, the image of T should cover the entire 
physical domain, that is, T should be onto. Also, the trans 
formation should be one to one. In terms of the grid, this 
means that the grid lines should not overlap the physical 
boundary and should intersect only at points corresponding 
to intersection points on the mesh in the computation domain 
Requiring T to be one to one and onto is equivalent 
to saying that the system T(s)=p must have one and only one 
solution in the computational domain for each point p in 
the physical domain. This provides the motivation for 
looking at the following general problem: 

Pick an open set DcR n , where R ' 1 is euclidean n-space, 
and let C be an open bounded set such that CcD. If 
F:DcR n ~R n is a continuous mapping and yeR n is given, how 
many solutions of F(x)=y exist in C? 

The difficulty in solving this problem lies in the 
fact that in general the solutions do not vary continuously 
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with F or y. This difficulty may be resolved by looking 
instead at the difference between the number of solutions 
for which the Jacobian of F is positive and the number of 
solutions for which the Jacobian of F is negative. Loosely, 
this is what is called the degree of F at y with respect 
to C. 


2 . 1 Defining the Degree of a Mapping 

A more precise definition of the degree of a mapping 
F takes on different forms depending on what restrictions 
are placed on F. What follows are essentially the defini- 
tions presented in references [S] and [0], 

2.1-1 Def i n i t i on . Let C / '° n be an open bounded set and let 
F:CcR n _ R n be continuously differentiable on C. Pick 
y^F(aC) and let r = lxeC|F(x) = y}. If F'(x) is nonsingular 
for all xer then one defines the degree of F'at y with 
respect to C by 

deg ( F , C , y ) = x | r sign det F'(x). 

In [0], Ortega and Rheinboldt actually define the 

degree in terms of an integral and then show that it has 

the equivalent form given above. 

On removing the restriction that det F'(x) t 0 for 

xer the definition becomes 

deg ( F , C , y ) = ! im deg( F, C , y. ) 
k— 

where lim y k = y and each element of { y ^ } 
k -»» 

satisfies y k { F ( a C ) and det F'(x) t 0 whenever F(x) = y k . 


Actually, one can make the stronger statement that 
for any such sequence {y k } there is a k Q such that 

deg(F,C,y) = deg(F,C,y k ) for k>k Q [0, p. 159]. 

The Weierstrass approximation theorem makes it 
possible to extend the definition of the degree of a mapping 
to a continuous function. 

2.1- 2 Definition . Let F:CcR n - R n be continuous on the 

bounded open set C. Define ||F||= | F(x)| where |«| 

is the Euclidean norm. Then for y^F(aC) one defines the 
degree of F at y with respect to C by 
deg(F.C.y) = 1 i m deg ( F . ,C,y) 

j— J 

where {F ^ } is a sequence of maps which are continuously 

differentiable on an open set D=>C and which satisfy 
lim || F i - F 1 1 = 0 

2 . 2 Properties of the Degree 

The principal properties cf the degree are given 
below. Excellent proofs may be found in [S], [0] and [ H ] . 

2.2- 1 Theorem . Let F:Ce. R n - R n continuous on the open 
bounded set C and letr = IxeC | F( x ) =y } . For any y^F(aC) 
there exists a quantity, deg(F,C,y), which has the propertie 
listed below. It i s : 

I • Integer va 1 ued 
2. Invariant under homotopy 
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If W:Cx[0, I ]cR n+1 - R n is continuous, then 
for any zeR n satisfying W(x,t) t z whenever 
(x,t)e aCx [0,1], deg ( W ( • ,t),C,z) is constant 
for all te [0, 1 ]. 

3. Dependent only on boundary values 

If 6:C«R n -R n is continuous and Gf 0C = F| qC , 
then deg(F,C,y)=deg(G,C,y). 

4. Invariant under translation 

For any zcR n . 

deg( F-z,C,y-z)=deg{ F, C,y ) . 

5. Invariant for points which can be connected by a 
continuous path avoiding F(aC) 

See Figure 1. 

6. Invariant unaer the excision from C of any closed 
set Q satisfying QAr = 0 

In other words, if QAr=0, then deg{F,C,y) = 
deg( F, C-Q,y) . In particular, if Q=C, deg( F, C-Q,y )=0. 
This property will be called the Excision Property. 
The Excision Property can be used to prove a very 
important result which is called the Kronecker Theorem in 

[0, p. 161]. 

2.2-2 Theorem (Kronecker). If F:CcR n - R n is continuous on 
the bounded open set C, yrfF(aC) and deg(F,C,y) ^ 0, then 
the equation F(x)=y has a solution in C. 

Proof : Suppose F has no solutions in C. Let Q=C. Since 

ytfF(Q), the Excision Property implies deg( F, C,y )=0. 


Q.E.D. 
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2.3 A Topological Definition of the Degree 

Dugundji [0] presents an alternate formulation for 
the degree of a mapping. He defines the degree of a 
mapping f : S — S where S is the unit n-sphere in R n , that is, 

S = {x c R n | | x | = 1 }. 

This degree can be shown to be equivalent to the analyti- 
cally defined degree in the previous sections. 

Before defining this degree t several terms must be 
discussed. 

2.3- 1 Definition . A set EeR n is called a linear variety 
if XpX 2 cE implies xx j + ( 1 -x ^eE for all real x. 

2. 3- 2 D efinition . A hyperplane in R n is an ( n- 1 ) dimen- 
sional linear variety. If n=l then a hyperplane will be a 
point. For n=2 it will be a line, and for n=3 it is a plane. 

2.3- 3 Definition . If {x Q ,Xp ,x n ) is a set of n+1 points 

in R n , then the convex hull is called an n-simplex . It will 
be denoted by 6 = (x Q ,Xj , . . . ,x n ). 

The points x Q ,Xj , . . . .x^ are called the vertices of the 
n-simplex. If the vertices lie on a hyperplane in R n , then 
the n-simplex is said to be degenerate . Now if (xj,...,xj) 

are the coordinates of point x^, then the volume of an 
n-simplex [F, p. 208] is given by 
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1 

n ! 

det(x 1 -x 0 , 

X2-Xq. • • 

* ’ x n” x o^ 



1 

n ! 


xLx 1 

xLx 1 

... x^x 1 



X 1 x o 

x 2 x o 

x n o 



det 

• 

• 

• 

• 

• 

• 




x n -x n 

x n -x n 

... x n -x n 




X 1 x o 

x 2 x o 

x n x o 



An n-simplex is degenerate if and only if 
det(x r x 0 ,x 2 -x 0 * n -V = 0. 

The next three definitions will be used to explain 
the term "ordered n-simplex." 

2. 3- 4 Def l n i t i on . A binary relation Aina set A is a 
subset acAxA, 

2.3- 5 Def i n i t l on . If a is a binary relation in a set A, 
then a is trichotomous if exactly one of the following is 
true for each x.yeA: 

x Ay , x=y , yAx. 

2.3- 6 Def i n i t i on . Let a be a binary relation in a set A. 
Then a is a total order if it is transitive and trichotomous 
[ G , p . 2 ] . 

2.3- 7 Def i n i t i on . An ordered n-simplex [D, p. 336] is an 
n-simplex together with a total ordering on its vertices. 

Therefore, if the vertices x ,Xj,...,x n of an n-simplex 
satisfy x Q <x j <. . . <x n , then "<" totally orders the set 
fx Q ,Xj , . . . ,x n J. Therefore, the n-simplex <5 = ( x Q , x j , . . . , y ) 
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is an ordered n-simplex. Such a simplex will be denoted 

[fi] * Cx Q ,Xj x n ]. The sign of the ordered simplex 

is the sign of det(xj-x .Xg-x t ...,x R -x ). 

Now suppose x 0 .Xj , . . . ,x R _ j is a set of n points on S 
having a diameter less than 1 so that the convex hull of 
the set does not contain the origin. Then the convex hull 
can be projected onto S by choosing the points on S lying 
on the directed rays which start at the origin and pass 
through the convex hull. The points on S form what will be 
called the spherical (n-1 ) - simplex 6 = ( x 0 * • • • » x n _| ) • 

The spherical simplex 6 is degenerate if and only if 
x ,Xj,...,x j, lie on a hyperplane in R n passing through 
the origin, that is, if and only if ( x Q , x j , . . . , x R _ j , 0 ) is a 
degenerate n-simplex in R n . An ordered spherical ( n- 1 ) - 
simplex is a spherical (n-1 )-simplex with a total order 
on its vertices. The sign of an ordered spherical 
(n-1 )-simplex [a] = [x ,...,x j] is defined to be the sign 
of the n-simplex [x Q x n _j,0] in R ° ^ D * p * 337 J- 

The next two definitions, which can be found in 
[D, p. 337], complete the terminology needed to define the 
Dugundji degree. 

2.3-8 Definition . A triangulation a of S is a decomposition 
of S into a finite number of nonoverlapping, nondegenerate 
spherical (n-1 )-simplexes such that each face of an (n-1 ) - 
simplex is the common face of exactly two (n-1 )-simplexes. 
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2.3- 9 Definition . Suppose S and e are unit n-spheres in 
R n . {Different symbols are used to make the concepts more 
clear.) Let a be a triangulation of S. A proper vertex 
map 91 a -z is a map defined only on the vertices of the 
spherical (n- 1 )-simplexes in a and is such that whenever 
x Q ,Xj , . . . ,x n _ j are vertices of a simplex in a, the set 

{9 ( x 0 ). <p( Xj ),..., <p ( x n _ | ) }cz has diameter less than 1. 

Under the proper vertex map tp : a— z there will be a 
unique simplex <p(o) lying on e corresponding to each simplex 
oea . There will be a unique ordered (n-1 )-simplex 
<p[o] = [<p(x 0 ),<p(x 1 ),...,<p(x n _ 1 )] on E corresponding to each 

ordered ( n - D- spherical simplex [o] . The sign of [o] may 

differ from that of <p[o]> and the family of sets {<p(a)|°ea} 
may not form a triangulation of e since it may contain 
overlapping simplexes and degenerate simplexes. However, 
the family does have the fundamental property presented in 
the following theorem which Dugundji proves [D, p. 237]. 

2.3- 10 Theorem . Suppose a is a tr i angu 1 ati on of S and 
<p:a -Ea proper vertex /nap. Let y be any point not on the 
boundary of any set <p(o). If p(y,a, 9 ) is the number of 
positive <p[<j] containing y ana n(y,a,<p) is the number of 
negative, then the number D(y, a, <p) = p(y, a, q>)-n(y , a, <p) is the 
same for all yeE not on the boundary of any <p(o). 

Since D(y, a, 9) is independent of y it can be denoted 


D( a, 9 ) . 
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Now if F:$- i is continuous then the compactness 
of S makes it possible to find a tri angulation a of S such 
that the diameter of F(o) is less than 1 for each oeA. 

Then if 9 p:A-z is the proper vertex map defined by 

9 P ( x ) = F ( x ) for each vertex x of a, Dugundji [D, p. 339] 

shows that the number D(A,9p), where 9p is the proper vertex 

map associated with a, is independent of the triangulation 
of S. He calls the quantity D(A, 9 p) the degree of F. 

Since a and 9 p actually depend only on F, D{ A, 9 p) can be 

denoted D(F). 

Like the analytically defined degree, this degree 
is invariant under homotopy [D, p. 239]. 

2.3-11 Theorem . If F:S - z i s homotopic to F:S^z, then 
D ( F ) = 0(F). 

Now let V be the unit n-ball in R n , that is, 

V * (x cR n | |x |<1}. Dugundji 's degree can be extended to a 
continuous map H:V-V provided H|^ maps S into S. S is 
clearly the boundary of V. Dugundji calls such maps regu 1 ar . 
The technique for determining the degree of H is analogous 
to what is done to obtain D(F) for F:S-S. V is triangulated 
into n-simplexes such that each face net on S is the face 
of exactly two n-simplexes. Then a regular vertex map is 
defined on the triangulation. 9 is a regular vertex map 
of a triangulation a of V if 9 maps each vertex on S to a 
point on S and 9 1$ is a proper vertex map. To calculate 
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the degree of H, which will be denoted D peg ( H ) , one chooses 
the regular vertex map defined by <p h (x) = H ( x ) for 

any vertex xca. Then choosing y e V-S such that y is not 
on the boundary of any q> H (o) one commutes 

D reg ( H ) = (number of positive <p H Co] containing y) 

- (number of negative q> H [o] containing y). 
D r eg ( H ) depends only on H, and if H and H are homotopic 
in such a way that the image of S remains on S throughout 
the entire deformation, then D reg (H) = D reg (H). Fur ther- 
more, Dugundji proves the following very useful result. 

2.3-12 Theorem . Suppose H : V — V is a regular map. Let 
F = H| S :S - S. Then D(F) = D reg (H). 

This theorem provides the information needed to show 
that Dugundji's degree is equivalent to the analytically 
defined degree. The following lemma will be used in the 
proof. 


2.3-13 Lemma . Suppose the following hypotheses are given: 

1. H:V-V is a continuously differentiable regular 
map 

2. yeV-S and r = (xeV | H(x ) = y } 

3. H '( x ) is nonsingular for all xer 

Then there exists a trUngulati on a of V with associated 
regular vertex map cp^ such that whenever oeA contains xer 
and is nondegenerate, 

sign [ o ] = sign det H'(x). 
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Proof : For each xcr choose a neighborhood N x of x so that 

either sign det H'(p)>0 for all points pcN x or sign 
det H '( p ) <0 for all points peN x . Choose the neighborhoods 
small enough so that the family of sets { N x | xc r > is dis- 
joint. Then triangulate V so that each N x contains a non- 
degenerate equilateral simplex.o x in which x lies, that 
is, a nondegenerate simplex in which the distance between 
any two vertices is the same. Call this triangulation a. 

Now suppose o y ei and c v = (x„,x, , . . . ,x_) with the vertices 
x x u i n 

labeled so that O x 3 = [x ,Xp...,x 3 is positive. Since 

<P H (x i ) = H(Xj), i = 0, 1 n, the sign of <p h I > x 3 

= sign det C H ( x j ) - H( x Q ) , . . . , H( x n ) - H(x Q )]. 

However, H(x i )-H(x Q ) = H'(x Q ) (x.-x Q ) + o(|x.-x 0 |) for 

i = l,...,n. Therefore, the sign of <p H [o ] = sign det 
[H '(x 0 )(xj-x 0 ) + o( |x j - x Q | ),..., H'(x 0 )(x n -x 0 ) + o(x n -x Q )]. 

Now if |x i - x Q | = e for i=o,l,...,n then expanding the deter- 
minant above yields det [H '(x ) (x^ -x , . . . ,x R -x ) ] plus terms 
of order o(e n ). Therefore, det [H '(x ) (Xj-x , . . . ,x n -x )], 
which has order 0(e n ), is the dominant term, and the other 
terms can be neglected. Consequently, the sign of 

<p h C<j x 3 = sign det [H' (x Q ) (Xj-x Q x n -x Q )3 

= sign [(det H'(x Q ))*(det (Xj-x 0 x x Q ) ] 

= sign det H ' ( x ) since [o ] is positive. However, 

v X 

since x Q eN x , sign det H ' ( x Q ) =• sign det H'(x). 


Q.E.D. 
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2.3-14 Theorem . Suppose H:V-V is a continuously differen- 
tiable regular map. Pick y e V-S and let r * { x c V | H ( x ) = y}. 
Let F =s H : S — S . If H'(x) is nonsingular for all x c r then 

D(F) > deg ( H, V, y ) . 

Proof : Since 2.3-12 says D(F) * ^reg-^’ ** suffices to 

show that D fe g(H) = deg (H,V,y). 

Choose a triangulation a of V as specified in 2.3-13 
and let be the regular vertex map associated with a. 
Without loss of generality, one can assume that q> H (o) is 
nondegenerate for all oea because any degenerate 9 H (o) can 
be approximated by a nondegenerate simplex 9 ( 0 ) where 9 is 
defined on all vertices p in a so that | <p ( p )-9 H ( P ) | < e 
for a given e. According to Dugundji [D, p. 338], D ( a , 9 ) = 
D ( A , 9 U ) = D ( H ) if e is sufficiently small. 

Furthermore, one can also assume that y does not lie 
on the boundary of any 9 H (o) for oea since property number 
5 of Section 2.2 implies deg(H.V.y) = deg(H,V,p) for all 
pcV-S. 

Therefore, it follows from 2.3-13 that D re g ( H ) = 
number of positive 9 ^( 0 ) containing y 

- number of negative 9 h (<j) containing y 

= 1 sign det h(x). 

X e r 

Q. E.D. 

The following corollary shows that the restriction 
that H'(x) be nonsingular for all xcr can be removed. 
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2.3-15 Corollary . Suppose H:V-V is a continuously 
differentiable regular map. Pick yeV-S and let 
r « (xeV |H(x)«y 1. Let F=H | s : S - S. Then D( F ) *deg( H. V, y ) . 
Proof : By Ortega and Rheinboldt [0, p. 159], there exists 

a sequence {y k 1 which converges to y and has the following 
properties : 

1 . Each y k *H( S ) 

2. For each y. , det H'(x) i 0 for all x such that 

r\ 

H ( x ) * y k 

3. For some k Q , deg( H, V , y ) = deg( H , V , y k ) for all 

k>k 
— o 

So pick k* so that k*>k Q and whenever k>k*, y k eV-S. Then 
by 2.3-14 D(F)=deg(H,V,y k *)=deg(H,V,y). 

Q. E.D. 

It should be noted that this corollary still holds 
if H maps some of the interior points outside of V. The 
points outside of V can be projected onto S so that one 
obtains a mapping from V into V. 

2.4 Applications to Grid Generation 

The usefulness of degree theory in grid generation 
surfaces when one studies a grid generating transformati on 
T. One might immediately note from the Kronecker theorem 
that determining the degree at every point in the physical 
domain would show whether or not T were onto. Unfortunately, 
the degree is not always easy to compute in practice. 
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One therefore looks Instead at how the degree can be used 
to prove some things about those quantities, such as the 
Jacobian of T, which can be easily computed. 

Recall that if AcR n , BcR n , then A homeomorphic to B 
means there is a continuous one to one, onto mapping from 
A to B whose inverse is also continuous. It is clear that 
T should be a homeomorphi sm from the computational domain 
onto the physical domain. 

The following result shows that if T is a homeomor- 
phism, its Jacobian does not change sign. 

In all of the theorems which follow l n =[0,l] n , JT = 
Jacobian of T, C°=interior of C, C c =comp 1 ement of C,aC = 
boundary of C and ncR n is homeomorphic to I p . 

2.4-1 Theorem . I , ‘ is a homeomorph i sm from I n to n and T 
is continuously differentiable, then the Jacobian, JT, of T 
has one sign in 1°, i.e., either JT(x)>0 for all xel° or 
JT ( x ) < 0 for all xel°. 

Proof : Suppose by way of contradi cti on that JT(x Q )>0 while 

J T ( x j ) <0 for some x 0 -XjeI°. Let y Q =T(x o ) and y 1 = T(x 1 ). 

Define p : [0, 1 ] - R n by 

p(t) - T ( (l-t)x 0 +tXj). 

Then p(0)=y Q , p(l)=y 1 , and p(tWT(aI R ) for t e [0,l]. Hence, 
by property 5, 1 =deg ( T, I R , y Q ) = deg ( T , 1° , y j ) = -1. There- 
fore, either JT(x)>0 or J T ( x ) < 0 for all * e I°. 


Q. E.D. 
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In an algebraic grid generation algorithm, the con- 
struction of T will be based on boundary information. 

The next theorem shows that requiring T to be a homemor- 
phism from the boundary of the computat ional domain to 
the boundary of the physical domain will insure that the 
image of T covers all of the physical domain. 

2.4- 2 Theorem . If T : I n — R n is continuously differentiable 

and T maps aI R homeomorphical ly onto an. then T(I )3 q . 
Proof: Let S be the unit n-sphere in R n . By Dugundji 

[D, p. 353], the Dugundji degree D of a map which is a 
homeomorph i sni from S to S is +1 or -1. But ai n and an 
are homeomor^ni c to S. Therefore, from 2.3-15 it follows 
that for any yen 0 , deg(T,I°,y) = ± 1. Therefore, by the 
Kronecker Theorem (2.2-2) n lies in the image of T. 

Q.E.D. 

Smith and Sritharan [SS] show that i f an additional 
hypothesis is added, one can obtain a much stronger conclu- 
sion: 

2.4- 3 Theorem . If T:I n -R n is continuously differentiable, 

T maps aI R homeomorphical ly onto an and JT(x) ^ 0 for all 

x el 0 , then T is a homeomorphi sm from ! n to n. 

The next theorem shows that the Jacobian changes sign 
when the image of T overlaps the physical boundary. Theorem 

2.4-3 and Theorem 2.4-4 show that it is important that T be 
constructed so that its Jacobian does not change sign. 
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2.4-4 Theorem . Suppose T : I n — R n has the following properties: 

1. T is continuously differentiable 

2. T maps al n homeomorphical ly onto an 

3. m( T( I n )— q)> 0 
Then JT has a sign change. 

Proof : Let C ( I n ) = (xeI n |JT(x) = 01 Sard's Theorem [0, 

p. 130] says that m( T( C( I n ) ) ) = 0. Since m( T( I p )-n)> 0, there 
exists z*e T ( I n )— n such that T(x) = z* implies JT ( x ) t 0. 

Now, choose we[T(I n )] c . By property 5, 

deg (T,I°,z*) = deg (T,I°,w) = 0. Since deg (T,I°,z*) = 

1 sign JT ( x ) , the Jacobian values at all x satisfying 
{ x | T ( x ) = z * } 

T ( x ) = z* must cancel each other. 

Q. E. D. 

2 . 5 Additional Topological Questions 

Section 2.4 suggests other questions which should 
be asked. Can a continuously differentiable homeomorph i sm 
from al n to an always be extended to a continuously 
differentiable homeomorph i sm from I n to n? If not, under 
what conditions is such an extension possible? How can 
one guarantee that a mapping from I n to R n will be a 
d i f f eomorph i sm? 

The answers to these questions will provide valuable 
information for creating an algebraic grid generation 
mapping. Although this paper does not answer all of these 
questions, partial answers were presented in the previous 
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section. Also, the following example shows that contin- 
uously differentiable boundary homeomorphi sms cannot always 
be extended. 

2 

2.5-1 Example . Suppose Tr^^R is continuously differen- 
tiable and maps the boundary of the square homeomorphical ly 
onto the boundary of the nonconvex region o shown in figure 
2. Let p be the point indicated and ap = /as \ If 

w 

T,(p) = _aI( D ) and Tp ( p ) = _aT(p) then 
as p a n 

T(P+AP) = T ( P ) + T'(p)Up) + o ( | Ap | ) 

= T(p) + a^T^P) - An ^2 ( P ) + o { | a P | ) 

When | a P | is small, the terms of order o ( | a P | ) are 

negligible in size when compared to a^Tj(p) and AnT 2 (p)* 

Therefore, those terms may be neglected from the equation 

above. However, then it is clear that T ( p + a p ) must lie 

outside the boundary of n. Consequently, T cannot be a 

homeomorph i sm from I- to n. This is illustrated in figures 

3 and 4 which show the result of attempts to construct a 

tensor product spline transformation that maps the square 

onto n. In each case points overlap the boundary near 

the "V" shaped corner. 

The first grid was obtained by choosing the B-spline 
coefficients so thac the transf ormat i on approximated a 
transfinite bilinear interpolation mapping. This is dis- 
cussed in Chapter 4. The second grid was obtained by chang 
ing some of the coefficients in order to minimize a func- 
tional which is described in the next chapter. 


dv+d 








3. 


AN ALGEBRAIC GRID GENERATION MAPPING 


In this chapter an algebraic grid generation tech- 
nique which uses a transformation consisting of tensor 
product B-splines is discussed. In the first section, 
finite difference approximations to the transformed deriva- 
tives of a first order partial differential equation are 
examined. The effect of the size of the Jacobian on smooth 
ness and orthogonality is discussed, and its influence 
on local truncation error is examined. The next section 
defines the particular transformation of interest in this 
paper and discusses the properties of the building blocks 

for this transformation: kth order B-splines. The final 

% 

section discusses a functional which can be used to modify 
the transformation so that the grid lines are distributed 
more smoothly and are nearly orthogonal at points of inter- 
section. 


3. 1 A First Order Example 


If C and n are the computational coordinates 
fying 0<e<l and 0<n<l, and x and y are the physica 
dinates, then the grid on the physical domain will 
of coordinate lines produced by a mapping 


, satis- 
1 coor- 
consi st 
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If u t « F(x,y,u,u x ,u y ) is a first order partial differen- 
tial equation defined on the physical domain, then the 
chain rule yields (u u ) * (u x u y ) x J where J * fx^ x n l, the 

L y e y J 

Jacobian matrix for the transformation T. Hence 
* u x u y ) * (u 5 u n ) x J' 1 

* < U E u n>* pn - x nj^/jT 

where JT * |J| = x^y n - Xrtfc. It is clear that the partial 

differential equation can be transformed once the elements 

of J are computed. These elements may be approximated 

by differences when explicit formulas are not available. 

The transformed expressions for u and u show immediately 

a y 

that the grid must be structured so that JT t 0 at all 
mesh points U,n)- 

Once the partial differential equations are trans- 
formed, difference approximations can be written for u^ 

and u . Large truncation errors in the approximations 
n 

will affect the solution of the partial differential equa- 
tions. One can obtain an expression for the truncation 
error at mesh point ( C i , n ^ ) by doing a Taylor series expan- 
sion at (5j.ni). If u^ * u(5 j ,n j ) , then 


u i + i , j ' u u + * u ee + u eee !#i 3 * H0T 

“i-l.J * U IJ * U E*E + “EE i#! 2 - U EEE I^i 3 + H0T 
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where HOT * higher order terms. Subtracting these two 
equations and solving for ug yields 

U E - u U1.1 - “1-1-1 ’ “EEC + H0T 

2 At 

Similarly, 

% ■ "1.1*1 • u l.t-l - u nnn + H0T 

TTrt 

Therefore 

u x * L (yn6gu-yg6 n u) - _1_ (y n u CC g( A g) 2 -ygU nnn (A n ) 2 ) 
JT 6 JT 

+ ... 

where a^u and a n u are the central difference' approximations 
for and u , respectively. The truncation error is 

(y n u ?:cC (A C ) 2 ‘ y C u nnn (A n )2) + *** 

Now If r * (*), then 

JT ■ Vn - x n^E 

* <r c xr n ) . (0,0, 1) T 

■ |r ( | |r n | sin e 

where o is the angle of intersection of the grid lir-s 
at (c.n). Again, the importance of JT i 0 is evident, 
but one can also see why the grid lines should be as ortho- 
gonal as possible. The expression for JT implies that 
the truncation error is inversely proportional to sin e. 
However, according to Thompson, Warsi and Mastin [TWM, 
p. 82] a departure from orthogonality of up to 45° is 
usually tolerable. 
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3.2 8-splines 


The napping T discussed in this paper has the forn 


T( C. r) 



m n 

i*l j«l c ‘ij P ij (5,n) 


m n 

i* 1 j«l e ij B ij (5,n) 


0<C<1 


0<n<l 


where the 8^, 1*1, ...m; j*l,...,n are tensor products 


of B-splines and the coefficients B.j, i*i,...,m; 

j*!,...,n are real numbers. In this section, the terms 
B-spline, spline function and tensor product B-spline are 
defined, and some of the Important properties of these 
functions are discussed. 


3.2- 1 Defining B-splines 

The following definition is from A Practical Guide 
to Splines by Carl de Boor [de B, p. 108]. 

3. 2- 1-1 Definition . If t * I tj 1 is a nondecreasing sequence, 

then the i-th normalized B-spline of order k for knot se- 
quence t is defined by 

B i,k t (xJ * ^ 1 i+k” t i )^i t i+k^"* ) +~ 1 where XeR * 

The sequence t may be finite, infinite or biinfinite. 
The expression [tj , . . . ♦tj + j l J(* -x) 1 denotes the kth 

k- 1 

divided difference of {• -x)* or the leading coefficient 

▼ * 

of the polynomial of degree k which interpolates 
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k- 1 

at tj t i+k* The not# * ion (•-*)+ represents the trun- 

ks 1 

cated power function (t-x)* which is defined by 

(t-x)*"* for t>x 

| 0 for t< x 

The * indicates that the kth divided difference above should 

k 1 

be evaluated by holding x fixed and considering (t-x)*’ 1 

as a function of t only. Nevertheless, since Bj k t (x) 
changes as one chooses different values for x, it is clearly 
a function of x. 

The definition above differs slightly from the original 
definition given by Curry and Schoenberg. Their B-spline 

M i k t is relatecl t0 8 1 k t th€ e< 3 uation 

M i,k,t 3 ^ k/ ^ t i+k“ t i^ B i,k,t ^ de B ’ p * 109 ^* 

3.2.2 Properties of B-spllnes 

A kth order B-spline Bj k t is a piecewise polynomial 

of degree k-1 with breakpoints at t i ,...,t i+k . On each 

interval (tj,tj +1 ), B t k t is a polynomial of degree k-1 
or less. For convenience it will be assumed that Bj k t 

is continuous from the right at breakpoints. 

B-splines have many properties which make them 
convenient for applications involving computers. One impor- 
tant property is their small support. If x4[tj,t i+k 3, 
k 1 

then (t-x)* will be a polynomial of degree k-1 or less 
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on Hence Ct 1 .....t Uk ] (*-x)J‘ l - 0. Therefore, 

B i,k,t * x * * 0 for x * ^ t i ,t t+k^ 

This implies that the support of Bj k t can lie in 
at most k intervals of the form Ctj 3. Therefore, 

if {B } represents the sequence of B-splines of order k 

for the knot sequence t * <t A >, it follows that only the k 

B-splines Bj_ k+ j, Bj -k+ g B j can have su PP° rt ln an y 

given interval [tj,t 4+ j]. 

The next two results, which are proved in [de B, 
p. 110] and [de B, p. 130], respectively, shbw that B-splines 
form a partition of unity, i.e., the sequence <B ^ consists 

of nonnegative functions which sum up to 1. 


3. 2. 2-1 Theorem . If (B ^ > is the sequence of B-splines of 
order k for a nondecreasing sequence t ■ ( t i >, then 


i B, ( x ) 

i 1 


V e,(x) 

i «p-k+l 


for any xe(t p jt q ) where p and q are such that p-k+1 and 
q+k-l lie in the index set for t. 


3. 2. 2-2 Theorem . If Bj is the ith element of the sequence 

of B-splines of order k for a nondecreasing sequence 
t * (tj). then B | ( x ) >0 for tj<x<t 1+k . 

One can think of the "B M in B-splines as representing 
the word 'basis," for when the knot sequence t is chosen 
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appropriately, the kth order B-spliner for t fore a basis 
for the piecewise polynomial space P b - . P. _ is the 

* * * t v *»{|V 

notation used by de Boor [de B, p. 100] to represent the 
space of piecewise polynomials of degree k-l which have 
breakpoint sequence c and which satisfy smoothness conditions 
specified by v. If k * then the nonnegative 

sequence v • (v^}* gives the number of smoothness conditions 
at each Cj , 1 * 2,...,m. For example, if v | * 3 then any 
fcP k,t v must have dt l* ast 3 smoothness conditions at 
S M , that is, the function, its derivative and second deriva- 
tive must be continuous at K r The dimension of P k ^ v 

is km-5 \»i . 
i »2 1 

The following theorem of Curry and Schoenberg [de B,C] 
shows how the knot sequence t should be chosen so that 
the corresponding B-spline sequence forms a basis for P. r . 

3. 2. 2-3 Theorem (Curry and Schoenberg). 

Let c * (Cj be a strictly increasing sequence and 

m 

v * (V| >2 be a nonnegative integer sequence such that 

v,<k for all i. Set n • k+1 (k- Vl ) ■ km-1 v , and let 
l ~ 1*2 1 i«2 1 

a ■ L 

t * (tjl" be a nondecreasing sequence such that 

(1) t i<t 2 i**‘<t| l <Ci an <* ?m+l5 t n+li***i t n+k 
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(ii) for i*2,...,m, the number Cj occurs exactly 
k- v times in t. 

Then the sequence Bp...,B n of B-splines of order k for 
the knot sequence t is a basis for r viewed as 
functions on C ^ ^ ^ n -»- 1 ^ ‘ 

This theorem shows how the number of knots at a break- 
point translates into the amount of smoothness there. 

Since the number occurs exactly k- v } times in t and 
v. represents the number of smoothness conditions at Sj, 
the number of smoothness conditions at equals k minus 
the number of knots at . Hence if k=4 and ^ , 2<j<m, 
occurs exactly once in t then the piecewise polynomials 
generated by will satisfy three smoothness condi- 

tions at 5 j , i.e., the piecewise polynomials, their first 
derivative and their second derivative will be continuous 
at . 

3.2. 3 Spline Functions 

In early studies of splines, a spline function of 
order k was defined to be a piecewise polynomial of degree 
k-1 with k-2 continuous derivatives. However, in this 
paper the more general definition in [de B] is used. 

3. 2. 3-1 Definition . If t = {t t > is a nondecreasing sequence, 
then a spline function of order k with knot sequence t is 


any linear combination of the B-splines of order k for 
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the knot sequence t. If one denotes the collection of all 
such functions by S k t then 

s k,t *?i B i .k.t= =*1 real for *“ “• 

It is clear that when t has the form described in the Curry 
and Schoenberg theorem 3. 2. 2-3, $ k ^ y on Ct| c ,t n+ j]. 

The first derivative of a spline function B i k t 

can be found by using the differences between successive 
coefficients. The following result, proved in [ae B, p. 
138], shows that the derivative of a spline function of 
order k will be a spline function of order k-l. 


3. 2. 3-2 Theorem . Let * ^ B i k t be a kth order spline 
function constructed with B-splines Bj k ^ corresponding 
to a nondecreasing sequence t = { 1 4 > . Then the first deri- 
vative of B i k t is given by 


d(V B. 


dx 


i i i ,k, t 


) = t ( k - 1 ) - a; . | 

t i+k-l' t i 


The value of a spline function f = * o- B, . t at a 

j j j i s i * 

point x satisfying ^ i <x<t^ + j is a convex combination of 
the k coefficients a j+i _ k • * • • » a j • For t i <x<t i+l’ then 
f(x) * j Ia j B j,k,t {x) = j r s j. k + i “j B j,k,t (x) Wlth the B j,k,t 
satisfying j* B j k t^ s * an<i B k^ x ^~° for al1 j* 
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B-spline coefficients model the functions that they 
represent. In other words, the coefficients are 
approximately equal to the yalue of the function at certain 
points. This is illustrated in the next section. 

Carl de Boor [de B] proves the following result con- 
cerning the relationship between a spline function and 
its B-spline coefficients. The notation ||f|i[ a denotes 

max If (x) I. 

xcCa.b] 


3.2. 3-3 Theorem. Let ?o. B. . be a kth order spline 
i l i,k,t 

function constructed with B-splines Bj k t corresponding 

to a nondecreasing sequence t = (tj). Then there exists a 

positive constant D k , depending only on k, so that for all 
i. 


l i' 1 u k 


|?«. 
J J 


j.k 


, t il 


tt, 


i + 1 


t i+k-l 3 


3.2.4 Variation Diminishing Splines 


Given an f known to lie in P k ^ v one can write it 

in the form f = £ The Curry and Schoenberg Theorem 

i = 1 

(3. 2. 2-3) shows how one obtains the B-spline basis and 
the following lemma suggests how one might obtain the 
coefficients. Its proof may be found in [de B, p. 116], 
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3. 2. 4-1 Lemma (de Boor and Fix). Let Bj be the sequence 

of B-$plines of order k for a nondecreasing sequence 
t * |t j }. Let Xj be the linear functional defined for all 

f by H f . j" 1 (-) k * I ' r ^ k ‘ 1 * r) (t i )f tr) (t i ) where 
r*0 

*(t) * (t i+J -t) ••• (t i+jc -t) and is some arbi- 

trary point in the open interval (Tptj^). Then 
x. Bj * for all j. 

Hence, if f * J °i B i 4t follows That a k , l<k<n may 
i*l 

be found by computing x k f * x k (za s a k* By explicitly 

writing out the expression for x^f one can easily show 
[de B, p. 159] that aj = f(t.) + 0(|t|) if t. is any point 

in (tj,t i+k ) and |t| * m J x ^i + l" 1 !** How ^ ver » if 

* tf\ 1 < i <.n where tf * (t i + | ♦ ••• + T i+k _j) y ( k -l) Then 

aj = f(tf) + 0(|tj 2 ). Choosing ajsfftj) for l<i<r yields a 

shape preserving approximation called Schoenberg's variation 
dimini shing spline approximation [de B, p. 159]. So if 
t* * {tp j the variation diminishing spline approximation 

to f, vf, is defined by 

vf = ; m*) b,. 

i *1 

This spline reproduces polynomials of degree one, i.e., 
if f is a straight line then vf * f. For any f the number of 
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times the spline approximation crosses a given line will 
be less than or equal to the number of times f crosses the 
line. From this it follows that if f is nonnegative, then 
vf is nonnegative and if f is convex then vf is convex. How- 
ever, since v has these shape preserving properties, it is 
not d very high order approximation. In fact, if g is a 
function defined on [a,b] and g has m continuous derivatives 
for some m>2, then de Boor [de B, p. 161] states that 
I Is- vg | «Cg k |t | , where c g k is a constant depending 

on the order of the spline function k and the function g. 

No matter how large m is, no exponent larger than 2 can be 
put in the inequality. De Boor shows that it is possible to 
obtain other spline approximations which are more accurate, 
but variation diminishing splines are convenient for appli- 
cations such as computer-aided design and grid generation 
where shape preservation is important. 

3.2.5 Tensor Product B-splines 

3. 2. 5- 1 Definition . Let R be the set of real numbers. If 

V is a linear space of functions mapping some set X into R 

> 

and W is a linear space of functions mapping some set Y into 
R, then for each veV and weW the tensor product , v*w of v 
and w is defined by 

vtw(x.y) * v(x)w(y) for (x.y)eXxY. 

Furthermore, the set of all finite linear combinations of 
the form vaw for some veV and weW is called the tensor 
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produc t, VIW of V with W. 

A typical element u of VIW has the form 

u -? vw 

where a^eR, VjeV, w j eW for j»l»....n. 

If V and W are the linear spaces of spline functions 
S h s and S k t , respectively, then the elements of VIW are 

linear combinations of tensor product B-splines. A tensor 
product B-spline B^ is defined by B i j(x,y)«B i ^ hiS ( x ) B j tk tt (y) 

where B. h is the ith B-spline of order h for the knot 

sequence s = fs,} and B 4 ^ is the jth B-spline of order k 

for the knot sequence t = {t,}. An element u of VIW will 

J 

be called a tensor product spline and will have the form 


where a^eR for ^11 i.j. When h=k=4 u may also be called 
a bicubic spline. 

Many of the properties of tensor product B-splines 
follow trivially from B-spline properties. For example, 
the tensor product B-spline B . j will be positive on its 

support since both Bj h $ and Bj k t are positive on their 
support. Furthermore, the support of Bjj is small. Since 

B i,h,s (x) * 0 for xrfSj ,s 1+ hJ and B J,k.t ( y> = 0 for y«t r tj +k 
it is clear that B i j(x,y) = 0 if eitherxtfSpS^] or 
y*ftj , tj +k ]. Hence the support of B^ lies in the shaded 
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area shown in figure 5. 

Tensor product B-splines also form a partition of 

unity. It follows from 3. 2. 2-1 that ieB (x,y) * 

ij ij 

IB (x) IB (y) = 1 for any (x,y) e ( s , s ) x ( t , t ) where p 
i i j j M M 

and q are such that p-h+1 and q+h-1 lie in the index set 

for sequence s and r-k+1 and m+k-1 lie in the index set 

for sequence t. 

Partial derivatives of tensor product splines are 
easy to compute since they reduce to derivatives of spline 
functions. 


a (El a B ( x , y ) ) 
*x ij ij ij 


9 B ( x ) B (y) 

W ij ij i j 

SB (y) d (i o B (x)) 
j j 3x i ij i 


a (i „ B ( x.y) ) 
ij ij ij 


a sz B ( x ) B (y) 
*7 i j i j i j 

SB (x) d (S B (y) ) 
i i 37 j ij j 


3.3 A Smoothing Functional 

The mapping T described in this paper uses tensor 
product B-splines to map the unit square onto a physical 
domain of arbitrary shape. This section shows that 
choosing the coefficients of the tensor product B-splines 
so that they minimize a certain functional can improve 
the quality of the physical grid produced by T. This 
functional is described and conditions under which it will 
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have a minimum are examined. 


3.3.1 Characteristics of the Functional 


The coefficients of the mapping defined by 


T(c.n) 




m n 

i »j j*l a i j B ij (5,n) 


m n 

i«l i*l B ij B ij {c,n) 


0<5<1 


0<n<l 


can be divided into two groups: boundary coefficients and 

interior coefficients. T uses the boundary coefficients, 
j« 1 , . . . ,m and ^ | id* i s "l*...,n to map the boun* 

dary of the square onto the boundary of the physical domain. 
Hence, the flexibility of their values is limited. The 
rest of the coefficients, the interior coefficients, can 
be moved around in order to change the characteri sti cs 
of the physical grid. To produce orthogonality in the grid 
’ines and maximize the smoothness of the distribution of 
grid lines one can choose the interior coefficients to 
minimize the functional 


' k “‘((If) 2 *(t^) 2 ) dA + /i 2 w 2 C0ot,^ 


JTU.n) = Jacobian of T at U.n) 


ox U.n) ox U.n) 
ot on 


0£ U.n) jal U.n) 
o? On 


where 
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« OX (C ,o ) 0^ (t,n) - ay ( C . n ) ox (t.n). 
"SI Tn on 


Dot ( 5 , n) ■ aT ( c, n) • oT ( t» n) 
oc On 


s 

ax (5, n) 


ax ( t« nil 


at 

• 

an I 


ajr (t.n) 


M < C. n) 


at 


an 


« ax ( c. n) ox (t.n) + 0£ (t»n) e£ ( e. n) 

ac On 05 On 

and Wj( 5 ,n), w 2 (t,n) = weight functions evaluated at (t,n). 

After the minimization of f is completed, where Wj is large 

the variation of the Jacobian values at nearby points will 
be small. Hence, Wj can be used to decrease skewness in 

a grid. Where w 2 is large, Dot will be small causing the 

grid lines to approach orthogonality. 

To avoid the tedious differentiation and integration 
of tensor product B-splines, the following discrete approxi 
mation to F can be implemented in computer algorithms: 

6 ’ L ( ( jt ^.j;V 2 * ( jt ..j»i; jt *j ) 2 W. 

J ' («)* (In ) 2 > 

+ P Q 2 

I I w,( Dot, , ) atari 

i = 1 .M 2 U 

where 

5 1<^ 2 < * * ’ <^p s 1 • 


0 


8 
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0 * <n2 < * * ,<f *q * 1 * 

JTjj * JT UjHj). Dotjj * D Qt (CfHj). 

ac « 1 / ( p- 1 ), An * 1 / ( q- l ), and 

the parameters w, and w 2 are weight functions. Both F 
and G depend only on the coefficients of the tensor product 
B-splines which compose T. 

Now 


m n 


Ox 

n 

(C.n) 

* E 
i * 1 
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J-l 

•u 
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OC 

(B U 

(C.n)) 
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n 





Ox 

7n 

(C.n) 

» E 
i * 1 
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j«l 

*u 
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an 

(e ij 

(C.n)) 
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n 





Hi 

(C.n) 
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(B, 1 
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i * 1 

J-l 

l J 

ac 

1 J 
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n 





Hi 

(C.n) 

* E 
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(B H 

(C.n)). 

tr\ 


i *1 

J-l 

1 J 

«n 

1 J 



Thus, for C,n fixed, JT(c.n) is a linear function in each 
coefficient i*l,...,m, J*l,...,n and Dot (c.n) 

is a quadratic polynomial in each coefficient. Since the 
terms involving Oot (C,n) and JT (e.n) are squared in G, 
one can see that G is actually a quartic polynomial in 
each coefficient. This suggests an elementary iteration 
method for finding the minimum of G: the cyclic coordinate 

method [B, p. 271}. 
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The cyclic coordinate method is a multidimensional 
search technique for minimizing a function of several vari- 
ables without using derivatives. It searches for a minimum 
along each coordinate direction. This method, when applied 
to a differentiable function, converges to a point where 
the gradient is zero [B, p. 273]. It can be applied to 
G if one treats each coefficient i"l m, 

j - 1 n as a variable representing a particular coordinate 

direction. This technique is discussed further in the 
next chapter. 

The importance of requiring that the Jacobian of T be 
of one sign was illustrated in Chapter 2. For this reason, 
if possible, the feasible region for the minimization 
problem is chosen to be a region where the Jacobian of T is 
nonnegative. Now since B-splines have small support, any 
given coefficient a fS or B fS only affects the Jacobian 
of T at a small number of points on the unit square mesh. 

By solving the inequality JT>0 for « rs at each of these 

points one can determine on what interval a fS must lie so 

that the Jacobian values at the points it affects are non- 
negative. This inequality is easy to solve since JT is 
linear in a fS . Repeating this procedure for each coefficient 

will, in most cases, produce a satisfactory approximation 
to the desired feasible region. However, since the boun- 
dary coefficients are fixed, there may sometimes be problems 
near the boundary. This is the case with the nonconvex 
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region examined in Section 2.5. The Jacobian will remain 
negative at one of its corner points even after the domain 
for the coefficients is restricted by using the procedure 
above. This is because the boundary points are fixed and 
not affected by the procedure. The Jacobian will also 
remain negative near this corner because of continuity. 


3.3.2 Convergence of the Smoothing Functional 


Under what conditions will the discrete smoothing 
functional 6 converge to a minimum value? Is it important 
that G be restricted to a region where the Jacobian of T is 
nonnegative? What happens if one of the tensor product 
coefficients becomes large? 

These are some of the questions which might be asked 
about G. The notation defined below will be used to discuss 
these problems: 

(A r >is a sequence in which each term represents a set 

2 

of coefficients for the mapping T : I ^ R defined by 
m n 

z £ **1 1 ® i i • n ) 

i = 1 j«l 


T ( C * n ) * 


t 0<6U. 0<n< 1 . 


m 

z 

i*l 


n 


z 

J-I 


«jj®ij(c*n) 


Each A f can be considered a discrete function defined 


by 
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A r (s.i J) 


•i r 

r 




1J 


if $ » 1 
If s t 2 


| i *1 | < i ( j * 1 


n. 


T r denotes the mapping obtained when the coefficients 

given by A f are used for T. 

JT denotes the Jacobian of T... 
r r 

IVU * 


It follows that if the sequence { A r > of coefficients 

converges to a single point then the corresponding values 
of 6 also converge. Hence, it is important to determine 
conditions which guarantee the convergence of the coeffi- 
cient sequence. Well, since the elements of {A r > can be 

viewed as points in R^ mn , the sequence converges if and 
only if it is a Cauchy sequence; however, a necessary condi 
tion for the convergence of { A r ) is that the sequence be 
bounded. The following theorem and corollary show how the 
Jacobian affects the boundedness of the sequence. 


3. 3. 2-1 Theorem . Suppose for all r T f maps el 2 homeomor- 
phically onto ea. If JT r (c.n)20 for all r and all points 
U.n) c I 2 , then either { A r > is bounded or JT r (c 0 .n 0 ) « 0 
for some point (c ,n 0 ) * I 2 . 
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Proof : By way of contradiction, suppose {A r } is not bounded 
and JT r U.n)>0 for all r and all points U.n) c I 2 . For 
any integer N there exists an A e (A } such that |A r L. V >N. 

i * I hi a a 

c. r N 

But this implies that either |>N or |p n |>N for some 

ij ij 

i,j. Now since n is bounded there exists M>o such that 

| p | <M for all pen. From3.2.3-3 it follows that for large 

enough N, max |T ( C • n ) | > 2M- 
0<e<i r N 
0<n<l 

Hence T r maps some point (fynj) e I 2 outside an. Now 

since JT„ >0 on l2, m(T (I 9 )-n )>0. But then 2.4.4 says 
r N c r N c 


that JT has a sign change. 
r N 


The corollary below follows immediately. 


Q.E.D. 


3. 3. 2-2 Corollary . Suppose for all r T f maps I 2 homeomor- 

phically onto an. If JT r (c.n)>0 for all r and all points 

U.n) e I 2 , then {A f > is bounded. 

One would like to show that the requirement JT r U,n)>0 

for all r and all points U.n) e I 2 is sufficient to guarantee 
the boundedness of { A r } . As indicated in 3. 3. 2-1, it is 
clear that if the magnitude of a coefficient is large enough, 
then the mapping T p associated with the coefficient will 
map some point in I 2 outside an. However, it is no longer 
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clear that m(T r ( I 2 )-q)>0 since JT r (g,n) may be 0 outside oo. 

Thus 2.4.4 cannot be used to obtain the contradiction that 
JT r must have a sign change as was done in Theorem 3. 3. 2-1. 
Although the writer has been unable to devise an acceptable 
proof to date, further study may show that the inequality, 
ro(T r ( l2)~o) > 0, is actually true. 


4. PROGRAM TENTEST 


This chapter discusses the computer program TENTEST 
which algebraically generates grids using tensor product 
cubic B-splines. A listing of TENTEST is given in the 
appendix at the end of this paper. 

The first section of this chapter presents the major 
steps involved in the computer algorithm. Sections 2 
through 5 examine the important features of the program, 
briefly discussing the subroutines involved. 

4.1 The Algorithm 

Although TENTEST contains almost a thousand lines 


of code, it is based on the following eight step algorithm: 


i. Input knot sequences {Sj } and ft i > consisting 
of values from [0,1]. 3 

ii. Compute the tensor product cubic B-splines 
corresponding to the knot sequences. 

iii. Choose initial coefficients to form a bicubic 
spline mapping from the square to a physical 
domain. 

iv. Use the mapping to plot a grid on the physical 
domain. 

v. If grid satisfactory, stop. If grid unsatisfac- 
tory, continue. 

vi. Input weights for smoothing functional. 

vii. Complete one iteration of minimization routine 
to obtain new coefficients. 
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viii. Go to step iv. 

There also exists a batch version of TENTEST which 
allows the user to request several iterations of the minimi- 
zation routine at a time. All the information needed to 
plot the initial and final grids is stored in files which 
can be interactively accessed after the execution of the 
program is completed. 

The programs were run on a PRIME 750 computer. The 
PRIMOS operating system, coupled with a PLOT 10 graphics 
package, was used to interactively draw the grids on a 
Tektronics 4014 terminal. The PRIME 750 can communicate 
at a baud rate of up to 9600 thus making it satisfactory 
for interactive graphics. 

4.2 Computing the Tensor Product B-splines 

Since B-splines are determined by the knots with 
which they are associated, the first concern of the user 
is to choose appropriate knot sequences. The user must 
pick two sequences s = { s A > and t= {t j > , placing them in 
file TENS0RDAT. The user actually picks only the "interior" 
knots for each sequence. In other words, he constructs 
two increasing sequences of numbers between 0 and 1. After 
reading the numbers from file TENS0RDAT, TENTEST places 
four 0's at the beginning of each sequence and four l's 
at the end of each sequence. By 3. 2. 2-3 {Curry and Schoenberg) 
and 3. 2. 3-1, the cubic B-splines associated with s and t 
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form bases for spline spaces S 4 s and S 4 t . The functions 

in each of these spaces will have three continuity condi- 
tions at each interior knot. The products of the 8-splines 
will form a basis for the tensor product of S, and S. + . 

The tensor product B-splines can be used to construct 
a transformation T on the square which maps the boundary 
of the square onto the boundary of a physical domain 
as described in Chapter 3. The user may obtain a better 
approximation to the boundary of the physical domain by 
increasing the number of interior knots in s and t or by 
redistributing the knots. This is discussed. in more detail 
in Section 4.5. 

On a given pxq mesh on the square with mesh points 
(£ u ,n y ). u=l,...,p, v=I,...,q, the values of the tensor 
product B-splines which compose T are fixed. Since these 
tensor product B-splines are the products of B-splines 
Bj,i=l, — ,m and Bj,j=l,...,n for some m and n, it is con- 
venient to store the function values and first derivatives 
of these B-splines at each and n y . Subroutine COMSPLINE 

uses the de Boor routine BSPLVD [de B, p. 288] to compute 
these values. BSPLVD calculates the function value and 
derivatives of all the nonvanishing B-splines at a given 
point. COMSPLINE stores the function values and first 
derivatives in two arrays: XSPLINE and YSPLINE. Therefore, 

after a call to COMSPLINE is completed, XSPLINE will con- 
tain the function value and first derivative of each B- 
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spline in iBi^.l at Cu« u«l,...,p and YSPLINE contains 
the function value and first derivative of each B-$pline 
in ffijljUl at n v ,v*l ,. . . ,q. Computing T or its partial 
derivatives at a mesh point becomes a matter of calculating 
the sum of the products of the tensor product coefficients 
with the appropriate elements of XSPLINE and YSPLINE. 

This computation is done in subroutine TENVALF. 

The next section explains how the coefficients are 
chosen initially. 

4.3 Choosing the Initial Coefficients 

Many different methods can be used to choose the 
coefficients initially. Since B-spline coefficients model 
the function they represent, one might simply choose the 
boundary coefficients to equal points along the boundary 
of the physical domain, and choose the interior coefficients 
to equal points known to lie in the interior of the physical 
domain. However, this creates the problem of deciding 
which interior points should be chosen as coefficients. 
Ideally, the original coefficients should produce a grid 
which is somewhat smooth so that only a few iterations 
are needed to obtain an acceptable degree of smoothness 
and orthogonality. 

For this reason, the computer program described in 

this paper initially selects coefficients which produce 

an approximation to the transfinite bilinear interpolant 

o 

of a mapping V : I 2 — R satisfying V:dl 2 - an. In reality 


3 



1 
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one need only define V on Mg* The user may provide para- 
metric equations which map the boundary of the square onto 
the boundary of the physical domain, or simply input a 
set of boundary points for the physical domain. In the 
first instance V is defined by using the parametric equa- 
tions. In the latter case V is obtained by linearly inter- 
polating between successive boundary points. The parametric 
equations below map the four sides of the unit square onto 
the four sides of the trapezoid as shown in figure 6. 


V(e.O) 
V(l.n) 
V( 5,1) 
V(0, n) 



The tra..sfinite bilinear interpolant U of V is defined 


by 

U( C, n) = ( 1- n)V( 5,0) + *V( 5,1) 

+ 5V(1, n) + d-5)V(0, n ) 

- (l-5)(l-n)V(0,0) - 5 ( 1 -n ) V ( 1 , 0 ) 

- (1-5) nV(O.l) - 5nV( 1,1). 

U agrees with V on the boundary of the square and hence 
interpolates V at an infinite number of points. Transfinite 
interpolants are discussed by William J. Gordon and Charles 
A. Hall in [G]. 

The program selects initial coefficients which pro- 
duce a variation diminishing spline approximation to U. 



Figure 6. Happing from computational domain to physical 
domain. 


5S 


Hence, if T is constructed from tensor products of B-splines 
B | « B | ^ s i * 1 ..... m and B j * B j ^ ^ j*l,...,n, which corres- 
pond to knot sequences s* { s ^ ^ j and t« {t j > j * 'respec- 

tively, then the initial coefficients of the tensor product 
splines are J^jj = U(s*,t*), i*l m; j»l n where 

s i*( s i + l , + * * * + s i + 3)/3* l = i m and tj = (tj^j + ...+tj^2)/3» 

j=l,...,n. Since variation diminishing splines yield exact 
approximations to linear polynomials, T will reproduce 
the boundary of any physical domain which can be divided 
into four line segments. Arbitrarily shaped boundaries 
can be approximated as accurately as desired by increasing 
the number of knots used to define the tensor product splines 
or by changing the placement of knots to increase the concen- 
tration in complex shaped areas of the boundary. 

The initial tensor product coefficients are constructed 
in subroutines BOUNCOEF and INNERCOEF. Figure 7 shows 
a grid on a trapezoid domain constructed with a mapping T 
having coefficients as described above. The grid is the 
image of T over an equally spaced mesh on the square. 

4.4 Minimizing the Smoothing Functional 

In TENTEST, the cyclic coordinate method is used 
to find the minimum of the smoothing functional G described 
in Section 3.3. As the name suggests, this method attempts 
to find the minimum of a multivariable function by cyclicly 
searching in the direction of each coordinate axis. For 
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6 , the coordinate directions are represented by the tensor 
product coefficients gjj, i*l m; j«l,...,n. 

The user must first decide what size mesh should 
be used to obtain a grid with acceptable smoothness and 
orthogonality. 6 is a function of 2 mn coefficients, however, 
since the boundary coefficients are fixed only 2 (m- 2 ){n- 2 ) 
coefficients are free. Therefore, in general, the mesh 
used for the minimization technique should contain at least 
2 (m- 2 )(n- 2 )points. 

The user must also decide on the size of the weights 
wj,w 2 for G. One can choose constant weights for both 
JT and Dot, or choose a weight function for Oot which pro- 
duces more orthogonality near the boundary of the grid 
than in the interior. Small constant weights of values 
between 1 and 10 can be used initially to determine how 
they affect the smoothness and orthogonality of the grid. 

Changing coefficient (or changes the value 

of the mapping T only on the support of the tensor product 
B-spline Bjj. Therefore, in order to locate the minimum 
of G in the direction represented by one need only 
consider the sum over those terms in G which contain the 
value of JT or Oot at mesh points ( 5, n) lying on the support 
of B . j . Subroutine CORANGE determines the range of summa- 
tion associated with each tensor product coefficient for 
a given mesh on the square, and function GF computes the 
sum over the range indicated by CORANGE. Figure 8 shows 
the support of a tensor product B-spline associated with 



Figure 8. Support for tensor Product 8-spline s 
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knot sequences s« (Sj and t« i t ^ > j - The shaded sec- 
tion represents the support of tensor product B-spllne 
Bg g. In order to minimize in the direction of coefficient 

5 it would be sufficient to look at the sum 


GF 


6 

r 

i *3 


7 
£ 

J *3 


V JT uij- J V 


i *4 j «2 




+ 67 9 

i 2 w 2 (Dot n raun. 
i *4 j*3 J 

Like G, the partial sum, GF, will be a quartic polynomial 
in each coefficient. 

All of this information is used by the minimization 
routine FFHIN. Each call to FFMIN produces one complete 
iteration of the cyclic coordinate method. For each coeffi- 
cient, the routine first determines the interval on which 
the coefficient must lie if JT is to be nonnegative at most of 
the mesh points affected by the coefficient. Then it calls 
eitner TESTMIN0, TESTMINL, TESTMINR, or TESTMINB depending 
on whether the interval is biinfinite, has a left endpoint, 
a right endpoint, or two endpoints. The chosen subroutine 
finds the location of the minimum of GF on the interval 
and changes the value of the appropriate coefficient accor- 
dingly. 


*•5 Distribution Functions 

If solutions of partial differential equations on 
a domain are to be accurate, the grid on the domain must 
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be concentrated in areas of rapid change such as boundary 
layers and shocks. In most cases concentration near the 
boundary of the domain can be easily accomplished through 
the use of distribution functions. 

Rearranging the points on the square mesh changes 
the distribution of grid points on the physical domain. A 
nonuniform distribution of points on the square mesh can 
be viewed as the Image of functions ?pl| - Ij. and 

93 : 1 j - ij defined on c and n« respectively. The grid is 


then generated by the mapping T defined by 


T(t,n) * TOf(c.n) 
where 1 - I 2 satisfies 


e (C»n) 


f j (€ ) 
fgU) • 


This is graphically illustrated in figure 9 . The grid on 
the physical domain is the image under f of an equally 
spaced mesh on the square. 

In the current version of TENTEST, the user may 
request one of tnree distributions for e and uniform, 
exponential, cr arctangent. Selecting the uniform option 
produces an equally spaced distribution. The distribution 
function is simply the identity function on I.. If the 
exponential option is selected, TENTEST calls routine 
EXPONENTIAL which maps cd. into *(c) « e cg -l 




63 


where c is a nonzero constant. If c>0, * concentrates 
the grid lines near the line corresponding to c*0. If 
c<0, ♦ concentrates the grid lines closer to the line 
corresponding to c*l. The grid in figure 10a was produced 
with tj(C) * 6 and f 2 (n) = *(n). The constant c is 4. In 

figure 10b, fj(6) * ♦(?) with c=5 and ^(n) = n. The 

degree of concentration increases or decreases as |c | is 
increased or decreased. In figure 10c, *j<C) = ♦(€) with 

c=2 and * 2 ( n) * n. 

TENTEST calls ARCTANGENT when the user selects the 
arctangent distribution option. ARCTANGENT maps c c I j into 

y(c) = arctangent (2cc-c) - arctangent { - c ) 
arctangent (c) • arctangent ( -c ) 

where c is a positive constant. This function concentrates 
grid lines near points corresponding to ?=0 and c = l simul- 
taneously. This is shown in figure lOd with <p j ( C) = S, 

9 ^( n) = r( n) and c = 5. 

Future improvements to TENTEST might include the 
addition of more distribution functions and the creation 
of a routine which allows the user to create his own dis- 
tribution function by interactively digitizing points on 
the unit square. The routine would then create a variation 
diminishing spline approximation to the points to form 
the distribution function. 

Since the distribution functions described in this 
section are defined on Ij, they can also be used to 



Figure 10c. Exponential dis- Figure lOd. Arctangent 

tribution on distribution 

5 with c = -2 on n with c*5. 




Figure 10. Concentrating grid points on trapezoid domain. 
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redistribute the knots which define the tensor product 
B-splines that form T. This will permit the user to con- 
centrate more knots in areas mapped to complex portions 
of the physical boundary so that T produces a better boun- 
dary approximation. Presently the user can choose to keep 
the original distribution on the knots or choose to redis- 
tribute the knots to obtain an exponential or arctangent 
di stribution. 


5. RESULTS AND DISCUSSION 


This chapter examines some of the grids produced by 
TENTEST. Physical domains of various shapes are illustrated. 
Some of the grids are for actual objects, such as an airfoil 
or part of the space shuttle, but most are simply grids on 
domains of various shapes and sizes chosen to illustrate 
the range of the program. 

The user's chief concern is the creation of an accept- 
able grid on a given physical domain in the shortest amount 
of time possible. Since the grid will be the image of a 
continuous mapping on the square, the best technique is to 
minimize the smoothing functional by using a grid generated 
from a coarse mesh. Then, once the new coefficients are 
obtained, the user can request that the grid be plotted 
using a much finer mesh. This technique is illustrated in 
the examples which follow. Most of the examples contain at 
least four grids: The image under T, with its initial 

coefficients, of a coarse square mesh; the image of a finer 
mesh; the image of the coarse mesh after several iterations 
of the minimization procedure; and the image of a finer mesh 
after application of the minimization procedure. Any other 
grids shown are chosen to illustrate grid concentration or 
other points of interest. In all the examples shown, only 
constant weight functions were used in the smoothing functional. 
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The first four examples show grids on domains with 
common geometric shapes: a trapezoid, a quadrilateral 

with unequal, nonparallel sides, a triangle and a circle. 
Since the domains are simply connected and convex, only a 
few interior points are needed for the sequences s and t 
which determine the tensor product B-splines that compose T. 

The next three examples show grids on domains which 
are not convex. The major concern with such grids is the 
overlapping of grid lines near the boundary. 

The last examples deal with grids around concr-te 
objects such as an airfoil or part of the space shuttle. 

The irregular boundaries of some of these grids make it 
necessary to use more knots to define T. 

For convenience, the following notation is used in 
this chapter. 

= number of B-splines Bj in the sequence corres- 
ponding to knot sequence s, or 4 + number of 
interior knots in s. 

N n = number of B-splines B, in the sequence corres- 
ponding to knot sequence t, or 4 + number of 
interior knots in t. 

wj s constant weight multiplied times the terms in 

the smoothing functional involving the Jacobian, 
JT, of T. 

wd ■ constant weight multiplied times the terms in 
the smoothing functional containing Dot. 
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Ng x N n will be the dimension of the tensor product 

spline space generated by Bjj * Bj x Bj , i«l N c ; 

j - I....N5. cpu = central processing unit - main control 
section of a computer. 


5. 1 Convex Domains 

The first three examples, which have linear boundaries, 
require only one interior knot for each of the knot sequences 
s and t. The simplicity of the domains also means that 
a very coarse grid can be used to minimize the smoothing 
functional. Four or five iterations produce good results. 

The circular grids in the fourth example require more 
interior knots. 

5.1.1 Trapezoid 

In this example Ng = N n = 5 and wj = wd = 1. The 
first picture in figure 11 is the grid obtained using the 
initial coefficients in Section 4.3. It is the image under 
T of an equally spaced 5x5 mesh on the square. This is 
the grid on which the minimization procedure was applied. 

Note that the number of grid points is 25, while the number 
of free coefficients is given by 2(Ng-2)(N n -2) = 18. 

Figure lib is a finer grid constructed using the same 
coefficients. Figure 11c shows how the initial 5x5 grid 
changes after five iterations of the minimization procedure. 
The new coefficients produce grid lines that appear to 
be nearly orthogonal at most grid points. The image under 
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Figure 11a Original grid. 



Figure lib Original grid 
refined. 




Figure 11c Optimized grid. Figure lid Optimized grid 

refined. 


Figure 11 Grids on trapezoid domain. 
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the new T of a 20x20 mesh is given in figure lid. The 
amount of cpu time used in the optimization process was 
1 minute and 25 seconds. 

In figure 12 the weights wj and wd have been changed 
to show what effect they have in the minimization process. 
Figure 12a shows how the initial 5x5 grid is changed after 
only three iterations when wj=0 and wd*l. Orthogonality 
is more pronounced, but the grid spacing is no longer as 
smooth. In the refined grid in 12b the spacing is very 
skewed near the top boundary. Figure 12c shows the 5x5 
grid after five iterations with wj=l and wd=0. The spacing 
is smoother but the grid lines are not orthogonal. Figure 
1 2d shows a finer grid. 

5.1.2 Quadrilateral with Unequal Sides 

Again, in this example = N n = 5 which means 
sequences s and t each contain one interior knot. Also, 
wj=wd=l. The minimization procedure was applied on the 
5x5 grid shown in figure 13a. Five iterations of the 
technique produced the grid in 13c. Figures 13b and 13d 
show refined versions of the grids in 13a and 13b, respec- 
tively. The five iterations of the minimization procedure 
required 2 minutes and 16 seconds of cpu time. In figure 
14 the optimized grids are concentrated near different 
parts of the boundary. In 14a an exponential distribution 
with parameter c*4 has been put on n. Figure 14b shows 
an exponential distribution on 6 and n with c*4 in each 
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Figure 12a wj«o, wd«l 



Figure 12b wj«o, wd«l 



Figure 12c wj»l, wd*o 
Figure 12 Effect of weights, 



Figur 1 2d wj«l, wd«o 


wj and wd 


Figure 13a Original grid 


Figure 13b Original grid 
ref i ned 




Figure 13d Optimized grid 
refined 


Figure 13 Grids on quadri lateral with unequal sides. 
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case. In figures 14c and I4d, an arctangent distribution 
with c-5 has been placed on £ and n, respectively. 

5. 1 . 3 Triangle 

In the previous examples, it was clear that each 
side of the unit square should be mapped to a side of the 
four-sided physical domain, but in the case of a triangle, 
which has three sides, this cannot be done. The boundary 
must be divided into four sections. The simplest thing 
to do is to divide one of the sides of the triangle Into 
two parts so that two sides of the unit square are mapped 
onto one side of the triangle as shown in figure 15. Figure 
16a shows the initial 5x5 grid constructed with * Nn * 5 
and wj«wd«l. Figure 16b shows a 20x20 grid constructed 
using the same coefficients. After five iterations of 
the minimization procedure, the initial 5x5 grid is trans- 
formed into figure 16c. Figure 16d shows a finer grid. 
Optimization required 2 minutes and 22 seconds of cpu time. 

5.1.4 Circle 

Variation diminishing splines reproduce straight 
lines exactly, but the same cannot be said about their 
approximation of nonlinear curves. For such curves the 
accuracy of the approximation depends on the number of 
knots used to define the spline function. For this reason 
more knots are needed to obtain a satisfactory mapping of 
the unit square onto a circular physical domain. For the 








77 


grids shown in figure 17, N$ * N n * 9 and wj*wd*l. Hence, 
there are five interior knots in both sequence s and 
sequence t. 

Note that 2( Ng-2) ( N n -2) = 98. Although this number 
indicates that a mesh of at least 98 points should be used 
for the minimization routine, the 8x8 grid shown in figure 
17a seems to produce an acceptable grid. One reason for 
this might be that the initial grid in 17a already appears 
to be quite smooth and orthogonal at most points. The 
major problems with orthogonality occur near the areas 
to which the corners of the square are mapped. These areas 
are indicated by the arrows in '7a. Figure 18a shows how 
the initial grid is changed after fifteen iterations of 
the minimization procedure. Figures 17b and 18b show finer 
grids. The fifteen iterations of the minimization procedure 
required 20 minutes and 14 seconds of cpu time. 

5.2 Nonconvex Domains 

The grids in this section show some of the difficulties 
in creating grids on domains which are not convex sets. 

5.2. 1 Nonconvex Quadrilateral 

Figure 19 shows the shape of the domain. This example 
was first mentioned in Section 2.5. The boundary of the 
unit square is mapped onto the boundary of the domain as 
indicated in figure 2. Example 2.5-1 shows that T will 
not map the square homeomorphical ly onto the domain even 




Figure 18a Optimized grid 



Figure 18b Optimized grid refined 

Figure 18 Grids on circular domain after optimization 
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after the coefficients are changed. This fact is supported 
by the negative Jacobian present at one of the corners of 
the square. The negative sign suggests that points near 
that corner will be mapped outside of the physical domain. 
This is confirmed by the grids illustrated. In this 
example, Nc = Nn * 5 and wj*wd*l. The 5x5 initial grid 
shown in figure 20a was used for the minimization procedure. 
The enlarged picture in figure 20b shows a finer grid. 

Figure 21a shows the result of four iterations of the 
minimization procedure. The nonnegative Jacobian require- 
ment pulls the grid lines into the interior of the domain. 
However, figure 22 shows an enlarged version of the corner 
which indicates that part of the grid still overlaps the 
boundary. This means that the minimization routine was 
unable to restrict all of the coefficients to intervals 
where the Jacobian of T is nonnegative. 

This is further indicated in figure 21b which shows 
a finer version of the grid in figure 21a. The four itera- 
tions of the minimization procedure required 1 minute and 
1 second of cpu time. 

5.2.2 Puzzle Pieces 

The next two domains, illustrated in figure 23, look 
like pieces from a puzzle. In each case N? * 19, Nn = 5, 
wj = 1 and wd=10. 

Grids on the first domain are shown in figures 24 
and 25. The minimization procedure was performed on the 



Figure 21a Optimized grid 



Figure 21b Optimized grid refined 

Figure 21 Optimized grids on nonconvex quadrilateral domain. 
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22x8 grid In figure 24a. Figure 24b shows a finer grid. 
Figure 24c shows the grid obtained after forty iterations 
and figure 24d shows a finer grid. The grids in figure 
25 show how the initial grid changes after two, five and 
fifteen iterations. The grid obtained after forty itera- 
tions is shown again for comparison. On this domain Tentest 
is able to pull all of the grid lines into the interior 
of the domain. 

Grids on the second domain are shown in figure 26. 

The initial 22x8 grid is shown in figure 26a and figure 
26b shows a finer grid. After forty iterations, the initial 
grid is transformed into 26c and a finer grid is shown in 
26d. Figure 27a shows a grid on the first domain concen- 
trated near the bottom boundary by using an exponential 
distribution on n with c=4. Figure 27b shows a grid on 
the second domain concentrated near the top by using an 
exponential distribution on n with c*-4. 

The forty iterations used for the first domain 
required 1 hour, 42 minutes and 23 seconds of cpu time, but 
the second domain required 2 hours, 4 minutes and 46 seconds 
for forty iterations. 

5.3 Grids for Specific Objects 

This section deals with grids about particular objects 
such as an airfoil. Tne boundaries often have peculiarities 
which make it difficult to obtain satisfactory grids. In 
many cases it may be difficult to maintain smoothness in 
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Figure 26a Original grid 


Figure 26b 


Original grid 
refined 





Figure 26c Optimized grid 


Figure 26d Optimized grid 
refined 



Figure 26 Grids on second puzzle shaped domain. 
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the grid while increasing orthogonal ity. Often the user 
must try to find an acceptable balance. He must also attempt 
to concentrate the grids in areas where rapid changes are 
likely to occur when partial differential equations are 
solved on the domain. 

5.3. 1 Airfoi 1 

The grids in this example are for the Karman - 
Trefftz airfoil. The parameters N$ = 19, N n = 9, wj*l and 
wd*.5. Hence, there are 15 knots in the s sequence and 
5 knots in the t sequence. Figure 28 shows how the domain 
can be viewed as having a boundary consisting of four parts. 
The minimization procedure was performed on the 21x12 grid 
in figure 29a. The grid lines appear to be orthogonal 
everywhere except near boundaries 1, 2 and 4. Note the 
sharp corners behind the airfoil. After one iteration the 
corners have been eliminated and the angles of the lines 
near the airfoil are not as acute. This is shown in 
figure 29b and in the finer grid in figure 30a. 

Solutions on a grid about an airfoil are usually more 
accurate if a higher concentration of points is placed 
near the airfoil boundary since this is the area most 
affected as air moves over the airfoil. Figure 30b shows 
a 30x30 grid concentrated near the airfoil boundary by 
using an exponential distribution on n with constant c*4. 

The minimization procedure required 6 minutes and 


36 seconds of cpu time. 
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Figure 29a Original grid 


Figure 29b Optimized grid 



Figure 29 Srids for KaVma'n - Trefftz airfoil. 
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Figure 30a Optimized grid refined 



Figure 30b Optimized grid concentrated near airfoil boundary. 
Figure 30 Optimized grids for Ka'rman - Trefftz airfoil. 



5.3.2 Spike-Nosed Body 
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According to [Sm, p. 130], the spike-nosed configura- 
tion occurs frequently in supersonic flow. R. E. Smith 
states that supersonic flow about such bodies is unsteady, 
with separation occurring near the nose-shoulder region. 
Therefore, the grids must be concentrated in that area 
[Sm, p. 48]. The boundary data for the grids shown in 
this section can be found in [Sm, p. 60]. Rotating the 
bottom boundary around a horizontal axis of symmetry produces 
a clearer picture of the actual body. The ratio of the 
length of the nose to the height of the shoulder is 2.14. 

As in the previous example, * 19, N n * 9, wj*l 
and wd=0.5. The 21x12 initial grid in figure 31a was used 
for the minimization procedure. Two iterations produce a 
small amount of orthogonality near the bottom boundary as 
shown in figure 31b. Additional iterations produce an 
undesirable wiggle in the grid lines near the shoulder. 

Figure 32a shows a finer grid and figure 32b shows a grid 
concentrated near the bottom boundary by using an exponen- 
tial distribution on n with c=4. 

The two iterations of the minimization procedure 
required 13 minutes and 1 second of cpu time. 

5.3.3 Shuttle 

The grids in this example are for a model of the 
space shuttle. The optimized grids are the result of ten 
Iterations on the 32x12 grid shown in figure 33. Para- 
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Figure 32b Optimized grid concentrated near boundary of 
spike-nose body. 



Figure 32 Optimized grids for spike-nosed body. 
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utters Ng * 29, N n - 9, wj*wd*l. Ten iterations of the 
minimization procedure produce a small amount of orthogona- 
lity near the boundary of the shuttle as shown in figure 
34. Figure 35 shows an optimized 32x20 grid concentrated 
near the shuttle boundary using an exponential distribution 
on n with c*4. The ten iterations of the minimization 
procedure required 1 hour and 43 minutes of epu time. The 
32x12 grid in figure 33 is the largest grid on which the 
minimization procedure has been applied. 
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Figure 33 Original grid for shuttle. 
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f 

Figure 34 Optimized grid for shuttle. 
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Figure 35 Optimized grid concentrated near shuttle boundary. 


6. CONCLUSIONS 


This paptr has examined an effective algebraic 
method for creating boundary fitted coordinate systems. 

The method, which involves a mapping T composed of tensor 
product B-splines allows one to regulate grid characteris- 
tics by adjusting the coefficients of the splines. Modifying 
the coefficients so that they minimize a smoothing functional 
enhances the smoothness and orthogonality of the grids 
generated by T. 

The method is implemented in the program TENTEST 
which gives the user control over the number and concentra- 
tion of grid points. T he user can also regulate the amount 
of smoothness and orthogonality in the grids by the selec- 
tion of weight functions for the smoothing functional. 

Suggestions for future revisions of TENTEST include 
the addition of more distribution functions to allow greater 
control over grid concentration. One might also investi- 
gate the possibility of adjusting the boundary coefficients 
during the optimization process so that the boundary points 
of the grid are affected by the minimization procedure. 

Ultimately, the true test for a grid comes when it 
is actually used to solve partial differential equations. 
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Therefore, the next stage of research must include solving 
problems on several grids produced by TENTEST. Then it 
•ay be possible to change the program into an adaptive 
technique which rearranges the grid points in response 
to gradient information from the evolving solution. 

Once these things are accomplished, one may attempt 
to use the technique to generate grids on more complicated 
multiconnected domains. This may involve the study of 
techniques for grid patching. 

Also, the Prime 750 computer is an excellent machine 
for graphics, but not very fast in solving problems involv- 
ing a large amount of computations. Hence, the possibility 
of creating a version of TENTEST which operates efficiently 
on a vector computer such as the VPS 32 at NASA Langley 
Research Center should be investigated. This will permit 
the user to run much larger and more complicated problems. 

Finally, once this grid generation technique has 
been thoroughly developed for two dimensional domains, a 
three dimensional technique can be attempted. T would 
become a mapping from the unit cube to the desired physical 
domain, composed of the tensor product of B-splines in 
the three coordinate directions. As in the two dimensional 
case, characteristics of the grid would be changed by 
changing the coefficients of the tensor product B-splines. 
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00002* PROGRAM TENTEST 

00003 tC 
, 00004 :c 

00005 :c TENTEST HAPS A SQUARE GRID <0.1>X<0,1> ONTO 

000061C A PHYSICAL DOMAIN OF ARBITRARY SHAPE THROUGH THE USE OF 
00007:C TENSOR PRODUCT B-SPLINES. THE ORIGINAL KNOT SEQUENCES 
00008 5 C MAY BE CHOSEN TO HAVE AN EQUALLY SPACED DISTRIBUTION* 

00009 JC EXPONENTIAL DISTRIBUTION* OR ARCTANGENT DISTRIBUTION* 
000101C SIMILAR CHOICES CAN BE MADE FOR THE DISTRIBUTION OF 
0001 1JC GRIDPOINTS ON THE SQUARE. 

000121C TENTEST CONSTRUCTS AN INITIAL GRID GENERATION MAPPING 

000131C CONSISTING OF A LINEAR COMBINATION OF TENSOR PRODUCT 
00014 SC B-SPLINES WITH THE COEFFICIENTS CHOSEN SO THAT THE MAPPING 
000151C YIELDS A VARIATION DIMINISHING SPLINE APPROXIMATION 
000161C TO THE TRANSFINITE BILINEAR INTERPOLANT OF A 
000175C FUNCTION WHICH HAPS THE BOUNDARY OF THE UNIT SQUARE 
00018:C ONTO THE BOUNDARY OF THE PHYSICAL DOMAIN. 

000191C IF THE USER REQUESTS A NEW GRID, TENTEST REARRANGES 

00020 1C THE COEFFICIENTS IN AN ATTEMPT TO MINIMIZE A FUNCTIONAL 
00021 1C G INVOLVING THE DIFFERENCE IN THE JACOBIAN OF THE GRID 
00022 1C GENERATION MAPPING AT ADJACENT MESH POINTS AND THE DOT 
00023 JC PRODUCT OF VECTORS TANGENT TO THE GRID LINES DN THE 
000241C PHYSICAL DOMAIN. 

00025:C 
00026 JC 

00027IC ROUTINES 
000281C 

00029 IC EXPONENTIAL 

00030 JC ARCTANGENT 

00031 JC FIXKNOTS 

00032 JC BOUNCOEF 

00033 JC INNERCOEF 

00034 1C COMSPLINE 

00035 1C TENVALF 

00036 tC TENSORVAL 

00037 JC JACOB 

00038 1C CORANGE 

00039 JC GF 

00040 1C FFMIN 

00041 IC CRIT 

00042 1C TESTMINO 

00043 1C TESTMINR 

00044 1C TESTMINL 

00045 1C TESTMINB 

00046 1C CUBIC 

00047 1C EXTREMES 

00048 1C NORM 

00049 1C 
00050 1C 
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00091 SC 
00052 JC 
00053 JC 
00054 :c 
00055 :c 
00054 JC 
00057 :c 
00058 I C 
0005? tC 

00040 :c 

00041 JC 
00042 JC 
00043 JC 
00044 JC 
00045 1C 
00044 !C 
00047JC 
00048SC 

00069 :c 
00070JC 
00071 :c 

00072JC 
00073* C 
00074 JC 
00075 JC 
00074 JC 
00077JC 
000785 C 
00079 JC 
00080 5 C 
00081 JC 
00082 JC 


THE FOLLOWING SUBROUTINES ARE ALSO REQUIRED* 

THEY HAY BE FOUND IN *A PRACTICAL GUIDE TO SPLINES' 

BY CARL DE BOOR, SPRINGER-VERLAG, 1978. 

BSPLVB... CONFUTES THE VALUE OF ALL POSSIBLE 

NONZERO B- SPLINES OF A GIVEN ORDER AT 
A GIVEN POINT* 

BSPLVD... COMPUTES THE VALUE AND DERIVATIVES 

OF ALL B-SPLINES WHICH DO NOT VANISH AT 
A GIVEN POINT 

INTERV.*. DETERMINES THE KNOT INTERVAL ON WHICH A 
GIVEN POINT LIES. ITS OUTPUT IS THE 
SUBSCRIPT IDENTIFYING THE KNOT WHICH IS 
IMMEDIATELY LEFT OF THE POINT. 

BVALUE * * . CALCULATES THE JDERIV-TH DERIVATIVE 

OF A SPLINE FUNCTION WHOSE COEFFICIENTS 
ARE STORED IN ARRAY BCOEF* THE VALUE OF 
JDERIV IS PROVIDED BY THE USER* 


TENTEST USES ROUTINES FROM A PL0T10 GRAPHICS 
PACKAGE TO PLOT THE GRIDS* 


VARIABLES 

NKNOTX,NKNOTY 

AND 


00083JC 
00084 JC 
00085 JC 
00084JC 
00087 JC 
00088JC 


NEWNOTX,NEWNOTY*. DIMENSIONS FOR SQUARE MESH* 

NX,NY.»» DIMENSION OF SPLINE SPACE IN X 

DIRECTION, Y DIRECTION. 

KX... QUANTITY OF NUMBERS TO BE ADDED TO THE FRONT 

AND BACK OF THE INTERIOR X KNOT SEQUENCE. 
ORDER OF B-SPLINES IN X DIRECTION* 


00089JC KY... 

00090JC 

00091JC 

00092 JC FNX,FNY, ♦ » 
00093 JC 

00094 JC BNX,BNY, , * 
00095JC 

00094 JC INX.INY* . » 
00097 JC 

00098 5 C INTX,INTY * • 
00099JC TX,TY.». 
00100 JC ALPHA, BETA... 
00101 JC 


QUANTITY OF NUMBERS TO BE ADDED TO THE FRONT 
AND BACK OF THE INTERIOR Y KNOT SEQUENCE* 
ORDER OF B-SPLINES IN Y DIRECTION. 

NUMBERS TO BE PLACED AT THE FRONT OF 
THE X AND Y KNOT SEQUENCES. RESPECTIVELY. 
NUMBERS TO BE PLACED AT THE BACK OF THE 
X AND Y KNOT SEQUENCES, RESPECTIVELY. 
DIMENSIONS OF INTERIOR X AND Y KNOT SEQUENCES, 
RESPECTIVELY. 

INTERIOR X KNOT SEQUENCE, INTERIOR Y KNOT SEQ. 

X KNOT SEQUENCE, Y KNOT SEQUENCE. 

ARRAYS CONTAINING COEFFICIENTS OF 
TENSOR PRODUCT SPLINE MAPPING. 
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00102 tC A*B 
00103 1C AND 

00104 1C AP* BP... ARRAYS CONTAINING COORDINATES FOR 

00105SC SQUARE MESH. 

00106 SC LEFTX, LEFTY... ARRAYS IDENTIFYING KNOT INTERVALS 

00107 SC ON WHICH SQUARE MESH COORDINATES LIE* 

00108 SC (LEFTX(I)*J IMPLIES TXC JK«A(IHTX< J+l) ) 

00109 SC W1*W2» . » WEIGHTS FOR JACOBIAN* DOT PRODUCT TO 

OOllOSC BE USED IN SMOOTHING FUNCTIONAL. 

00111SC Xf Y* ♦ ♦ TWO-DIMENSIONAL ARRAYS CONTAINING COORDINATES 

00112.C OF GRID POINTS TO BE PLOTTED. 

00U3SC KOUNTE... VARIABLE USED TO COUNT ITERATIONS OF 

001141C MINIMIZATION PROCEDURE* OR 

00115SC NUMBER OF CALLS TO ROUTINE FFHIN. 

00116SC 
001 17 SC 

001181C AUTHOR: BONITA VALERIE SAUNDERS 

00119SC DATE: JULY 1985 

00120 :c 

00121SC 
00122 SC 
001231C 

ooi24:c*m*t**mmm 

001 25 SC 

001261C 

00127 : COMMON/COEF/ALPHA (100,100), BETA (100,100) 

00128S COSWON/KNOTS/TX ( 100) , TY( 100) 

00129S COMMON/PARAM/FKOUNT 

00130 S C0MM0N/PARAM2/A(100)*B(100),NX,NY,KX,KY 

00131*. * , LEFTX (100), LEFTY (100) 

001325 CONMON/KNOT/NKNOTX,NKNOTY 

00133*. COMMON/WEIGHTS/W1 ,W2 

00134 *. COMHON/SPLINES/XSPLINE (50,100,2), YSPLINE (50,100,2) 

00135J COMMON/RANGE/IF IRST ( 100 ) , ILAST ( 100 ) , JF IRST ( 1 00 ) 

00136* * , JLAST (100) 

00137S REAL X(100,100),BCOEF(100),INTX(100),INTY(100) 

00138: t *Y ( 100,100) ,AP( 100) ,BP( 100) 

00139S CHARACTERS NAME 

00140.* INTEGER«2 STRINGS 28) 

0014 i: INTEGER42 NUMB, DATE (3) 

00142.* 1NTEGER42 TIME,TIHE1*TIME2 

00143: EQUIVALENCE ( STRINGS 1), DATE) 

00144*. EQUIVALENCE (STRING(4) ,T1HE) 

00145: EQUIVALENCE (STRING(5) , TIMED 

00146S EQUIVALENCE (STRINGS 7) ,TIME2) 

00147: NUMB 8 28 

00148! CALL TIMDAT(STRING,NUMB) 

00149S . WRITES 1,111) DATE 

00150: 111 FORMAT (3A2) 

00151 S WRITE(1*222) T1ME»TIME1,TIME2 

00152S 222 FORMAT (16, 16 *16) 
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00153: 


Pl-3* 14159 

00154: 


0PENd2,FILE-'TEN8ORDAT'> 

00155: 


0PEN!13,FILE-'NEWDATA') 

00154: 


OPEN! 14, FILE-' 0R0RID') 

00157: 


0PEN<16»FILE-'0RIS2') 

00158: 


Ul-0 

00159: 


82-0 

00140! 


KOUNTE-O 

00161 : 


PRINT*, 'INPUT NKN0TX,NKN0TY' 

00142: 


R£AD!l»t) NKNOTX,NKNOTY 

00163: 


NSAVEX-NKNOTX 

00144: 


NSAVEY-NKNOTY 

00165: 


PRINT*, 'WHAT IS KX' 

00144! 


READ(l,t) KX 

00147! 


PRINT*, 'WHAT IS KY' 

00148: 


READ< 1,«) KY 

00169: 


READd2,*> FNX,BNX»FNY,BNY 

ooi7o: 


READ<12,«) INX,INY 

ooi7i: 


R£AD<12»t) (INTX!I) ,I-1»INX) 

00172! 


READ! 12, t) (INTY(I) ,I«1,INY> 

00173! 


PRINT*, 'DISTRIBUTION FOR X KNOTS' 

00174! 


PRINT* ,' 1 -EQUALLY SPACED , 2-EXPONENT I AL , 3-ARCTANGENT ' 

00175! 


READ(1,«) KODEX 

00176! 


IF<K0DEX-2> 5,10,15 

00177! 

10 

CALL EXPQNENTIAL!INTX,INX»2. ) 

00178! 


GOTO 5 

00179! 

15 

CALL ARCTANGENT (INTX,INX, 5.) 

00180! 

5 

PRINT*, 'DISTRIBUTION FOR Y KNOTS' 

00181! 


PRINT*, '1-EQUALLY SPACED , 2 “EXPONENTIAL , 3 “ARCTANGENT ' 

00182! 


READd,*) KODEY 

00193! 


IF<K0DEY-2> 16,18,19 

00184! 

18 

CALL EXPONENTIAL! INTY,INY, 2.) 

00185! 


GO TO 16 

00184! 

19 

CALL ARCTANGENT! INTY,INY, 5.) 

00187! 

16 

CONTINUE 

00188! 


NX-INX+KX 

00189! 


NY-INY+KY 

00190! 


CALL F IXKNOTS ( KX ,KY , FNX , FNY , BNX , BNY , 

00191! 

* 

I NX , INY , INTX , INTY ) 

00192! 


CALL BOUNCOEF ( KX , KY , NX , NY > 

00193! 


CALL INNERCOEF < KX , KY , NX , NY ) 

00194! 


WRITE(16,*> !<ALPHA!I,J) ,J-t,NY),I-l,NX) 

00195! 


WRITE! 16,*) < <BETA(I, J>, J-l ,NY> ,1*1 ,NX) 

00196! 


DO SO I-l ,NKNOTX 

00197! 


A ( I ) -FLOAT < I - 1 ) / ! NKNOTX- 1 ) 

00198! 


DO 20 J-l ,NKNOTY 

00199! 


B < J ) -FLOAT < J- 1 )/< NKNOTY- 1 ) 

00200! 


IF(A(I) .GE.1.0) A!I)-. 99999 

00201! 


IF(B( J) «GE«1 .0) B!J)-. 99999 

00202! 

20 

CONTINUE 

00203! 

50 

CONTINUE 
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002041 


PRINT** 'DISTRIBUTION FOR COMPUTATIONAL X' 

00203; 


PRINT** ' 1 ■EQUALLY SPACED* 2-IXPONENT I AL, 3-ARCTANGENT' 

00204; 


READ(1*») KODEXX 

00207; 


IF(K0DEXX-2) 22*24*26 

00208! 

24 

CALL EXPONENT 1AL< A *NKNOTX >2.) 

00209; 


GOTO 22 

00210! 

26 

CALL ARCTANGENT ( A * NKNOTX * 5 . ) 

00211! 

22 

CONTINUE 

00212! 


PRINT** 'DISTRIBUTION FOR COMPUTATIONAL Y' 

00213! 


PR INT* , ' 1 -EQUALLY SPACED * 2-EXPONENT I AL , 3- ARCTANGENT ' 

00214! 


READ(1*S) KODEYY 

00213! 


IF(K0DEYY-2) 30*32*34 

00216! 

32 

CALL EXPONENT 1 AL < B * NKNOT Y * 2 . ) 

00217! 


GOTO 30 

00218! 

34 

CALL ARCTANGENT (B, NKNOT Y *5.) 

00219! 

30 

CALL COMSPLINE 

00220! 


CALL CORANGE (NKNOTX, NKNOTY) 

00221! 


DO 70 1-1 * NKNOTX 

00222! 


DO 60 J-l, NKNOTY 

00223! 


CALL TENVALF < ALPHA , LEFTX d ), LEFT Y < J ), KX , KY , I , J , 

00224! 

* 

X<I,J),0,0) 

00225! 


CALL TENVALF ( BETA *LEFTX ( I ) , LEFTY ( J > , KX , KY * I , J 

00226! 

* 

,Y(I*J) *0*0) 

00227! 

60 

CONTINUE 

00228! 

70 

CONTINUE 

00229! 


PRINT*, 'JACOBIAN YES(l) OR N0(0>' 

00230! 


READ(1**> JCODE 

00231! 


IF(JCODE.EG.l) GOTO 700 

00232! 


PRINT** 'COMPUTE DERIVATIVES YES(1> OR N0<0)' 

00233! 


READU,*) KCODE 

00234! 


IF(KCODE.EO.O) GOTO 80 

00235! 


PRINT*, 'INPUT DERIVATIVES DESIRED FOR X*Y DIRECT' 

00236! 


READ<1,*> JDX, JDY 

00237! 


PRINT*,' X Y X COMP DERIV Y COMP DE 

00238! 


DO 600 I 1-1, NKNOTX 

00239! 


DO 500 JJ*1, NKNOTY 

00240! 


CALL TENVALF ( ALPHA * LEFTX d I ) , LEFT Y < J J ) , KX , K7 * 1 1 , J J , 

00241 ! 

* 

XD, JDX, JDY) 

00242! 


CALL TENVALF < BETA , LEFTX ( 1 1 ) , LEFT Y ( J J ) , KX , KY , 1 1 , J J , 

00243! 

* 

YD, JDX, JDY) 

00244! 


PRINT*, A(II),B(JJ),XD, YD 

00245! 

500 

CONTINUE 

00246! 

600 

CONTINUE 

00247! 


GOTO 80 

00248! 

700 

CALL JACOB < NX, NY ,KX,KY,A,B) 

00249! 

80 

CONTINUE 

00250! 


CALL EXTREMES < X , Y , TMAX , TMIN , NKNOTX , NKNOTY > 

00251! 


CALL NORM < X , Y , TMAX , TM I N , NKNOTX , NKNO TY > 

00252! 


PAUSE 

00253! 


NN*2 

00254! 


NAME* 'PRODUCT GRID' 


ill 


002551 


WRITK13,*) NAME 

002541 


WRITE(13»*)NKN0TX,NKN0TY»NN,NN 

002571 


WRITE* 13, *> *X*I,1>, 1*1, NKNOTX) 

00258: 


WRITE* 13,0 (Y(I ,1) ,1*1 , NKNOTY) 

0025?: 


WRITE(13,0 *X*I, NKNOTY), 1*1, NKNOTX) 

00240: 


WRITE* 13,0 *Y*I,NKN0TY) ,I*1,NKN0TX) 

00241 : 


WRITE(13,0 X ( 1 , 1 ) , X ( 1 , NKNOTY ) 

00242: 


WRITE* 13,0 Yd, l),Y*l, NKNOTY) 

002431 


WR1TE*13,0 X ( HKNOTX , 1 ) , X ( NKNOTX , NKNOTY ) 

00244: 


WRITE* 13,*) Y < NKNOTX, 1),Y( NKNOTX, NKNOTY) 

002451 


CALL 1NITTC940) 

002441 


CALL TWIND0<0, 740, 0,740) 

002471 


CALL DWIND0*-,07,1.07,-,07,1.07> 

002481 

338 

DO 200 1*1, NKNOTX 

002491 


CALL MOVER (X(I,1),Y(I,1)) 

002701 


DO 100 J*l, NKNOTY 

002711 


CALL DRAWA*X(I,J),Y(I,J)) 

002721 

100 

CONTINUE 

002731 

200 

CONTINUE 

002741 


DO 400 J*l, NKNOTY 

002751 


CALL H0VEA(X(1,J) ,Y*1 , J) ) 

002741 


DO 300 1*1, NKNOTX 

002771 


CALL DRAWA(X(I, J) ,Y(I , J) ) 

002781 

300 

CONTINUE 

002791 

400 

CONTINUE 

002801 


WRITE* 14,*) NKNOTY, NKNOTX 

002811 


WRITE* 14, *) **X* I, J), 1*1, NKNOTX), J*l, NKNOTY) 

002821 


WRITE*14,*> (<Y(I,J), 1*1, NKNOTX), J*l, NKNOTY) 

002831 


CALL MOVABS* 0,760) 

002841 


CALL ANMODE 

002851 


PRINT*, ' ITERATION' ,KOUNTE 

002841 


PRINT*,' NX*', NX,' NY*', NY 

002871 


IF*K0UNTE»EQ»0) GOTO 410 

002881 


PRINT*, 'JACOBIAN WEIGHT*' ,W1 

002891 


PRINT*, 'ORTHOG WEIGHT*' ^W2 

002901 


PRINT*, 'OPTIMIZED 0N',N8AVEX,' BY' ,NSAVEY, ' 

002911 

410 

KOUNTE-KOUNTEtl 

002921 


PRINT*, 'DO YOU WANT TO CHANGE THE GRID, YES OR N0*0) 

002931 


READ(1,*) KODE 

002941 


IF*K0DE»EQ»0) GOTO 339 

002951 


PRINT*, 'CURRENT WEIGHTS ARE',W1,W2 

002941 


PRINT*, 'NEW WEIGHTS, YES(l) OR N0*0)' 

002971 


READ(1,*> KW 

00298: 


IF*KW»EQ«0) GOTO 401 

002991 


PRINT*, 'INPUT WEIGHTS FOR JACOB, ORTHOG' 

003001 


READ(1,*> W1,W2 

003011 

401 

CONTINUE 

003021 


NKNOTX*NSAVEX 

00303! 


NKNOTY-NSAVEY 

003041 


CALL FFMIN*ERMAX) 

003051 


PRINT*, 'OUTPUT JAC0BIAN,YES*1) OR N0*0>' 


GRID' 


00306! 


READd,*) JKODE 

00307: 


XF(JKOK.EQ.O) GOTO 3?? 

00308: 


CALL JACQB<NX,NY,KX,KY,A,B> 

00309: 

3?? 

CONTINUE 

00310 : 


PRINT*, 'CHANGE HUM OF MIDPOINTS, YES(l) OR N0<0>' 

oo3ii: 


READd,*) NODE 

00312: 


IF(KODE.EQ.O) GOTO 501 

00313: 


PRINT*, 'ENTER NUMBER OF MIDPOINTS FOR X DIRECTION 

00314: 

* 

Y DIRECTION' 

00315: 


READd,*) NEWNOTX , NEWNOTY 

00316: 


NKNOTX-NEWNOTX 

00317: 


NKNOTY-NEWNOTY 

00318: 


GOTO 502 

003,'?: 

501 

CONTINUE 

00320: 


NENNOTX-NKNOTX 

00321: 


NEWNOTY-NKNOTY 

00322: 

502 

CONTINUE 

00323: 


CALL ERASE 

00324! 


DO 450 I«l, NEWNOTX 

00325: 


DO 425 J-l»NEWNOTY 

00326: 


AP< I ) -FLOAT < 1-1 ) / ( NEUNOTX-1 ) 

00327: 


BP< J > -FLOAT < J-l )/< NEWNOTY-1 > 

00328! 


IF(APd)»GE,l*0) APU)». 99999 

0032?: 


IF(BP( J) .GE. 1.0) BP( J)-*99?9 

00330! 

425 

CONTINUE 

00331 : 

450 

CONTINUE 

00332! 


PRINT*, 'DISTRIBUTION FOR COMPUTATIONAL X' 

00333! 


PRINT* 1 -EQUALLY SPACED , 2-EXPONENTIAL , 3-ARCT AN ' 

00334! 


READd,*) K0DE3X 

00335: 


IF(K0DE3X-2) 452,454,456 

00336! 

454 

CALL EXPONENTIAL(AP, NEWNOTX, 2.) 

003375 


GOTO 452 

00338: 

456 

CALL ARCTANGENT <AP, NEWNOTX, 5.) 

00339! 

452 

CONTINUE 

00340! 


PRINT*, 'DISTRIBUTION FOR COMPUTATIONAL Y' 

00341: 


READ(1,*> K0DE3Y 

00342! 


IF(K0DE3Y-2) 458,460,462 

00343: 

460 

CALL EXPONENTIAL (BP, NEWNOTY ,2# ) 

00344! 


GOTO 458 

00345! 

462 

CALL ARCTANGENT ( BP , NEWNOT Y , 5 . ) 

00346! 

458 

CONTINUE 

00347! 


DO 480 >t, NEWNOTX 

00348: 


DO 470 J- 1, NEWNOT Y 

00349! 


CALL TENSORVAL ( ALPHA , NX , NY , KX , K Y , AP ( I ) , 

00350! 

* 

BP( J) ,X(I, J) ,0,0) 

00351! 


CALL TENSORVAL ( BETA , NX , NY , KX , KY , AP ( I ) , BP < J ) , 

00352! 

* 

Y(I, J), 0,0) 

00353! 

470 

CONTINUE 

00354! 

480 

CONTINUE 

00355! 


CALL EXTREMES CX , Y , TMAX , TMIN , NEWNOTX , NEWNOTY ) 

00356! 


CALL NORM ( X , Y , TMAX , TMIN, NEWNOTX , NEWNOTY ) 
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003571 PAUSE 

0035SI SOTO 331 

0035ft 339 CONTINUE 
00360: CALL FINITT<0,760> 

00361: PRINTS, 'KX IS ',KX 

00362: PRINTS, 'KY IS SKY 

00363: PRINTS, 'DISTRIBUTIONS' 

00364: PRINTS, ' 1 ■EQUALLY SPACED, 2 -EXPONENTIAL, 3- ARCTANGENT' 

003651 PRINTS, 'X KNOT DI8T, IS',KODCX 

00366: PRINTS, 'Y KNOT BIST. I8SK0DEY 

00367: PRINTS, 'COMPUTATIONAL X D1ST IS.' ,KODEXX 

00368: PRINTS, 'COHPUTATIONAL Y DIST IS',KODEYY 

00369: CLOSE( 12) 

00370 : close: 13) 

00371: CL0SEU4) 

00372: CALL TIHDAT < STRING, NUHB ) 

00373: WRITE<1,222) TIHE,TIME1,TIME2 

00374: STOP* 

00375: END 


00376:C 
00377 :c 
00378: 
00379 :c 
00380 :c 
00381 :c 
00382 :c 
00383 :c 
00384 :c 
00385JC 
00386:c 
00307 1 C 
00388 :c 
00389, C 
00390 :c 
00391 :c 
00392:c 
00393: C 
00394: 
00395: 
00396: 
00397J 
00398: 
00399: 
00400: 
0040i: 
00402 JC 
00403 :c 
00404 :c 
00405: 


SUBROUTINE EXPONENT I AL ( X , N , AC ) 

THIS ROUTINE PRODUCES AN EXPONENTIAL DISTRIBUTION OF 
POINTS BY SUBSTITUTING AN ORIGINAL SET OF NUMBERS 
U LYING BETWEEN 0 AND 1 INTO THE EXPONENTIAL 
FUNCTION (EXP(ASU)-l • >/<EXP(AC)*l) WHERE AC IS A 
PARAMETER WHOSE VALUE IS SUPPLIED BY THE USER, 

VARIABLES 

X, » .THIS AN ARRAY WHICH ON INPUT CONTAINS THE ORIGINAL 
SET OF NUMBERS AND ON OUTPUT CONTAINS THE EXPONENTIAL 
DISTRIBUTION OF NUMBERS. 

N...8IZE OF ARRAY X 

AC.. PARAMETER IN EXPONENTIAL FUNCTION 

REAL X(100) 

DO 10 1-1, N 
U-XCI) 

X(I)-(EXP(ACSU)-1.)/(EXP<AC)-1.) 

IF (X(I).GE.l.O) X<I)-. 99999 
10 CONTINUE 
RETURN 
END 


SUBROUTINE ARCTANGENT (X»N, AC) 
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004061C 

00407 l C THIS ROUTINE PRODUCES AH ARCTANGENT DISTRIBUTION 
00408 1C OF POINTS BY SUBSTITUTING AN ORIGINAL SET OF NUMBERS 
0040? I C U LYING BETWEEN 0 AND 1 INTO THE ARCTANGENT FUNCTION 
004101C <ATAN(AC>-ATAN(-ATAN<-AC)>/<ATAN(AC>-ATAN(-AC)) 

00411 1C WHERE AC IS A PARAMETER WHOSE VALUE IS 
004121C SUPPLIED BY THE USER. 

004131C 

00414SC VARIABLES 
004151C 

00416SC X. * .THIS AN ARRAY WHICH ON INPUT CONTAINS THE ORIGINAL SET 
00417IC OF NUMBERS AND ON OUTPUT C0NTAIN8 THE ARCTANGENT DISTRIBUT 

00418.C OF NUMBERS. 

004191C N» • .SIZE OF ARRAY X 

00420 1C AC.. PARAMETER IN ARCTANGENT FUNCTION 

00421tC 

00422 SC 

00423 S REAL X<100) 

00424* DO 10 I-1'N 

004251 U«X<I> 

00426 S X < I ) ■ < ATAN ( 2 . *AC*U-AC > -AT AN ( -AC ) ) 

00427J * /<ATAN(AC)-ATAN<-AC)> 

00428 S IF(X(I>.GE.1.0) X<I>*. 99999 

0042?S 10 CONTINUE 

004301 RETURN 

00431S END 

004321C 

00433 1C 

00434 1C 

004351 SUBROUTINE FIXKNOTS<KX,KY,FNX,FNY,BNX, BNY.INX. 

004361 * INY , INTX.1NTY) 

00437 1C 

00438 1C THIS ROUTINE PLACES KX COPIES OF FNX AT THE 

004391C BEGINNING OF THE INTERIOR X KNOT SEQUENCE. 

00440 1C KY COPIES OF FNY AT THE BEGINNING OF THE 
00441 1C INTERIOR Y KNOT SEQUENCE. KX COPIES OF BNX AT THE 
004421C END OF THE INTERIOR X KNOT SEQUENCE AND KY COPIES 
00443 1C OF BNY AT THE END OF THE INTERIOR Y KNOT SEQUENCE. 

00444 1C 
004451C 

00446 1C INPUT 

004471C 
00448 1C 
0044?!C KX... 

004501C 
00451 1C 
00452 SC KY... 

00453 1C 
00454 1C 

004551C FNX, FNY... 

00456 SC 


QUANTITY OF NUMBERS TO BE ADDED TO THE FRONT 
AND BACK OF THE INTERIOR X KNOT SEQUENCE. 

ORDER OF B-SPLINES IN X DIRECTION. 

QUANTITY OF NUMBERS TO BE ADDED TO THE FRONT 
AND BACK OF THE INTERIOR Y KNOT SEQUENCE. ORDER 
OF B-SPLINES IN Y DIRECTION. 

NUMBERS TO BE PLACED AT THE FRONT OF 
THE X AND Y KNOT SEQUENCES, RESPECTIVELY. 
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00437 1 C BNX,BNY... NUMBERS TO BE PLACED AT THE BACK OF THE 
00430 tC X AND Y KNOT SEQUENCES* RESPECTIVELY » 

00457 1C IMX*1NY.«» DIMENSIONS OF INTERIOR X AND Y KNOT SEQUENCES* 
00440 1C RESPECTIVELY* 

0044 tlC INTX»INTY* *« INTERIOR X KNOT SEQUENCE* INTERIOR Y KNCT SCO. 

00442 tC 

00443tC 

00444 !C OUTPUT < IN COMMON) 

00445 1C 

00464 1C TX»TY... X KNOT SEQUENCE* Y KNOT SEQUENCE 

004471 

00448 1C 

00447 1C 

004701 COMMON/KNOTS/TX (100 ) , TY ( 100 ) 

004711 REAL INTX(1)*INTY(1> 

004721 NX-INX+KX 

004731 NY-INYHCY 

004741 DO 100 I«KX+1*NX 

004731 J*IHCX 

004741 TX(I)«INTX< J) 

004771 100 CONTINUE 

004781 DO 200 I«KY+1,NY 

004771 J-I-KY 

004801 TY<I)«INTY<J> 

004811 200 CONTINUE 
004821 DO 5 I«1,KX 

004831 TX(I)«FNX 

004841 INDEX-HNX 

004831 TX< INDEX )*BNX 

004841 3 CONTINUE 

004671 DO 4 I-1,KY 

004881 TY<I)«FNY 

004871 INDEX-IfNY 

004701 TY (INDEX )«BNY 

004711 4 CONTINUE 

004721 RETURN 

004731 END 

00474 1C 
004731C 
00474 1C 

004771 SUBROUTINE BOUNCQEF<KX»KY*NX*NY> 

0047ft C 

00477 1C THIS ROUTINE COMPUTES THE BOUNDARY COEFFICIENTS 

00500 1C FOR TWO TENSOR PRODUCT B -SPLINE 8 <X,Y COMPONENTS). 

00501 1C SPECIFICALLY* IT COMPUTES ALPHA(I*1)*BETA(I,1) AND 

005021C ALPMA(I*MY) *BETA< I *NY) FOR 1*1 TO NX1 

00303 1C ALPHAd* J> *BETA(1*J) AND ALPHA(NX* J) *BETA(NX*J) 

00304 1C FOR >1 TO NY. 

00503 1C COEFFICIENTS ARE CHOSEN SO THAT THERE IS A 
00304 1C VARIATION DIMINISHING APPROXIMATION ALONG THE 
90307 1C BOUNDARY. 


00SO81C 

ooso?:c INPUT 
00510 !C 

0051 1!C KX,KY*.. ORDER OF SPLINE IN X DIRECTION, Y DIRECTION 
005121C NX»NY » . •DIHENSION OF SPLINE IN X DIRECT, Y DIRECT* 
00513 iC 

005141C DOTH COEFFICIENT SEQUENCES WILL HAVE DIHENSION 
00515: C NXtNY 
00516 :c 
00517*0 


005181C 


OUTPUT 

005195C 



00520 :c 

(IN CQHMON) 

0052i:c 

BOUNDARY COEFFICIENTS PLACED IN ALPHA, BETA ARRAYS 

00527 :c 



00523: 


CDHHON/COEF/ALPHA < 100 , 100 ), BETA < 100, 100 ) 

00524! 


C0HM0N/KN0TS/TX(100) ,TY(100> 

00525! 


DIHENSION TXSTAR(IOO) ,TYSTAR(100) 

00526! 


81X<T>=2.*T«. 

00527! 


G1Y(T)*0» 

00528! 


G2X(T)=3.+T 

0052?! 


G2Y(T)=2.*T 

00530! 


G3X(T)=4»tT 

00531! 


G3Y(T>=2* 

00532! 


G4X(T>*1 *-T 

00533! 


G4Y(T)=2.*T 

00534! 


PI*3. 14159 

00535! 


DO 100 1=1, NX 

00536! 


SUM=0 • 

00537! 


• DO 50 J=1,KX-1 

00538! 


SUH*SUH+TX(I+J) 

00539! 

50 

CONTINUE 

00540! 


TXST AR ( I ) =SUH/ ( KX - 1 ) 

00541! 

100 

CONTINUE 

00542! 


DO 200 J=1,NY 

00543! 


SUH=0. 

00544! 


DO 150 K=1,KY-1 

00545! 


SUH*SUH+TY ( J+K ) 

00546! 

150 

CONTINUE 

00547! 


TYSTAR(J)=SUH/(KY-1) 

00548! 

200 

CONTINUE 

00549! 


DO 300 1=1, NX 

00550! 


A*TXSTAR( I ) 

00551! 


ALPHA(I,1)=G1X(A) 

00552! 


BETA<I,1)=G1Y(A> 

00553! 


ALPHA ( I , NY ) *G3X ( A ) 

00554! 


BETA(I,NY)=G3Y(A) 

00555! 

300 

CONTINUE 

00556! 


DO 400 J=1,NY 

00557! 


B«TYSTAR(J) 

00558! 


ALPHA ( NX , J ) =G2X ( B ) 


005591 

00540: 

00541: 

00542: 

00543: 

00544: 

00545:c 

00344 :c 

00547:c 

00548: 

00549:c 

00570:c 

00571 :c 

00572 :c 

00573 :c 
00574 :c 
00575 SC 
00574SC 
00577 :c 
00578 :c 
00579:c 
00580 SC 
00581 SC 
00582 SC 
00583SC 
00584 SC 
00585 SC 
00584SC 
00587 SC 
00588SC 
00589 SC 
00590SC 
00591 SC 
00592 SC 
00593SC 
00594 SC 
00595SC 
00594 SC 
00597 S 
00598 S 
00599 S 
00400 S 
00401 S 
00402 S 
00403 S 
00404 S 
00405 S 
00404 S 
00407 : 
00408 S 
00409 : 


BETA(NX,J)«G2Y(B) 
ALPHAU,J>«84X<B> 
BETA(1,J)*G4Y(B> 
400 CONTINUE 
RETURN 
END 


SUBROUTINE INNERCOEF(KX,KY,NX,NY) 

THIS ROUTINE COMPUTES THE INNER COEFFICIENTS 
FOR TWO TENSOR PRODUCT B -SPUMES <X,Y COMPONENTS) 
SPECIFICALLY, IT COWHITES ALPHAS I, J) ,BETA<I,J) FOR 
1*2 TO NX-1, J*2 TO NY-1. 

THE COEFFICIENTS ARE COMPUTED THROUGH THE USE OF 
TRANSFINITE BILINEAR INTERPOLATION. THE 
INTERPOL ANTS ARE EVALUATED AT POINTS SO THAT THE 
RESULTING COEFFICIENTS PRODUCE A VARIATION 
DIMINISHING SPLINE APPROXIMATION TO THE 
TRANSFINITE BILINEAR INTERPOL ANT. 

INPUT 

KX,KY... ORDER OF B-SPLINES IN X DIRECTIONS DIRECTION 
NX, NY... DIMENSION OF SPLINE SPACE IN X DIRECTS DIRECT 

BOTH COEFFICIENT SEQUENCES HILL HAVE DIMENSION 
NX*NY 

TX»TY< IN COMMON)... KNOT SEQUENCE FOR X DIRECTS 
DIRECTION 


OUTPUT (IN COMMON) 

INTERIOR COEFFICIENTS PLACED IN ALPHA, BETA ARRAYS 


COMMON/COEF/ ALPHA ( 100 , 1 00 ) , BETA ( J 00 , 1 00 ) 
COMMON/KNOTS/TX < 100 > , TY ( 100 ) 

DIMENSION TXSTAR(IOO) »TYSTAR(100) 
GIX(T)«2.*T+1. 

G1Y(T)*0» 

G2X(T)*3»+T 

G2Y(T)»2.*T 

G3X(T)*4.»T 

G3Y(T)*2. 

G4X(T)*1.-T • 

G4Y(T)*2,tT 

FX<X,Y)*61X<X)*(1.-Y)+G3X(X)*Y 
« +G2X(Y)*X+(1.-X)*04X(Y) 
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*-GlX<0>t<l.-X>t(l.-Y>-62X(0>»X*(l#-Y> 
*-G3X ( 0 > t< 1 . -X > *Y -G2X (1 . > *XtY 
FY<X,Y)*G1Y<X)*<1.-Y>+G3Y<X>*Y 
$ +G2Y(Y)*XHt.-X>tG4Y<Y> 
t -G1Y<0)*<1.-X)*<1.-Y>-G2Y(0)*X*<1#-Y> 
*-G3Y<0>t<l.-X>*Y-G2Y<l#>«X*Y 
PI*3#14159 
DO 100 1*1, NX 
SUM*0. 

DO 50 J*i,KX-l 
SUM*SUM+TX(I+J) 

CONTINUE 

TXSTAR<I)*SUM/<KX-1> 

> CONTINUE 

DO 200 J*t»NY 
SUH*0. 

DO 150 K*1»KY-1 
SUH=SUM+TY<JFK> 

> CONTINUE 

TYSTAR( J)=SUM/(KY-1 ) 

) CONTINUE 

DO 400 1=2, NX-1 
DO 300 J=2, NY-1 
A*TXSTAfi(I> 

B=TYSTAR( J) 

ALPHA<I,J)=FX(A,B> 

BETA(I, J)*FY( A,B) 

) CONTINUE 
) CONTINUE 
RETURN 
END 


SUBROUTINE COMSPLINE 

THIS ROUTINE COHPUTES AND STORES THE FUNCTION 
VALUES AND FIRST DERIVATIVES OF ALL THE NON- 
VANISHING B-SPLINES AT EACH POINT OF A SQUARE 
MESH# 

IN ADDITION, IT DETERMINES THE KNOT INTERVAL 
ON WHICH EACH MESH COORDINATE LIES. 


INPUT 

(IN COMMON) 

A, 8... ARRAYS CONTAINING COORDINATES FOR SQUARE 

MESH, POINTS OF EVALUATION FOR B-SPLINES. 
NKNOTX,NKNOTY,.# NUMBER OV ELEMENTS IN A,B, 

KX,KY» « . ORDER OF B-SPLINES IN X DIRECTION, Y 

DIRECTION. 

TX, TY. # . X KNOT SEQUENCE FOR B-SPLINES, Y KNOT 
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00 * 41 :c 

00**2tC 

00*6310 

006*410 

006*510 

oo***:c 

00**7 1C 

oo**a :c 

00*6910 
00*70 :c 
00*7110 
00*7210 
00*7310 
00*7410 
00*7510 
00*7610 
00*7710 
00*7810 
00*7910 
00*8010 
00*8110 
00*8210 
00*8310 
00*8410 
00*8510 
0068*10 


SEQUENCE FOR B-SPL1NES. 

NX, NY... DIMENSION OF SPLINE SPACE IN X DIRECTION, 

Y DIRECTION 

OUTPUT 
(IN COHHON) 

XSPLINEtYSPLINE. .THREE DIMENSIONAL ARRAYS CONTAINING 

FUNCTION VALUES AND FIRST DERIVATIVES OF 
B-SPLINES IN X DIRECTION, Y DIRECTION AT 
EACH ELEMENT OF A,B* THE FIRST SUBSCRIPT 
IDENTIFIES THE B -SPLINE, THE SECOND 
SUBSCRIPT REPRESENTS THE POINT OF EVALUATION 
AND THE THIRD SUBSCRIPT (1 OR 2) INDICATES 
WHETHER THE VALUE REPRESENTS A FUNCTION 
EVALUATION OR DERIVATIVE EVALUATION. HENCE 
XSPLINE(3,2,1) WILL CONTAIN THE VALUE OF THE 
B-SPLINE IN THE X DIRECTION, B<3), AT A(2). 

LEFTX, LEFTY... ARRAYS IDENTIFYING KNOT INTERVALS ON 

WHICH MESH COORDINATES LIE. LEFTX<3)=4 WOULD 
MEAN THAT A<3> LIES BETWEEN TX(4) AND TX<4+1) 


REQUIRED ROUTINES! 
BSPLVD 
BSPLVB 


00687 1C 
00*8810 
00*89! 
00*90! 
00691! 
00*92! 
00*93! 
00*94! 
00*95! 
00*9*1 
00*97 1C 
00*9810 
00*9910 
00700! 
00701! 
007021 
00703! 
00704! 
00705! 
0070*1 
007071 
007081 
007091 
00710! 
00711! 


INTERV 

C0MM0N/PARAM2/A< 100) ,B< 100) ,NX,NY ,KX,KY »LEFTX( 100) 
* »LEFTY< 100) 

COMMON/SPLINES/XSPLINE<50,100,2),YSPL1NE(50,100,2) 
COMMON/KNOTS/TXUOO) ,TY(100) 

COMMON/KNOT /NKNOTX , NKNOTY 
REAL DBIATX(4,2),W0RK(4,4) 

IDER1M-2 

JDERIV*2 

INITIALIZATION 


1 

2 

3 


NUMX«NX+KX 
NUMY*NY+KY 
DO 3 1*1, NX 
DO 2 J*l, NKNOTX 
DO 1 KK=1,2 
XSPLINE < 1 , JfKK>=0. 
CONTINUE 
CONTINUE 
CONTINUE 
DO 9 1*1 , m 
DO 8 .1*1, NKNOTY 
DO 7 KK»1,2 




a -, i <£. 


00712* 

00713: 

00714: 

00715: 

00716: 

00717: 

00718: 

00719*. 

00720: 

00721 

00722 

00723 

00724 

00725 

00726 

00727 

00728 

00729 

00730 

00731 

00732 

00733 

00734 

00735 

00736 

00737 

00738 

00739 

00740 

00741 

00742 

00743 

00744 

00745 

00746 

00747 

00748 

00749 

00750 

00751 

00752 

00753 

00754 

00755 

00756 

00757 

00758 

00759 

00760 

00761 


YSPLINE(I,J,KK>*0. 

7 CONTINUE 

8 CONTINUE 

9 CONTINUE 

00 25 I*1,NKN0TX 

CALL INTERV(TX»NUMX,A(I) »LEFTX(I) ,NFLAB> 

CALL BSPLVDt TX , KX , At I ) , LEFTX ( I ) .WORK, OBI ATX , I DERI V ) 

00 252 JJ*1,2 

00 251 11*1, KX 

IB*LEFTX(I)-KX«I 

IF(IB.LE.O) GOTO 251 

XSPLINE(IB,I,JJ)*DBIATX(XI,JJ) 

251 CONTINUE 

252 CONTINUE 
25 CONTINUE 

PRINT*, 'DC ALL X SEO EQUAL ALL Y SEG YES(1> OR N0(0) 

REAO(l,t) KODE 

IF(KODE.EQ.O) GOTO 28 

DO 27 JJ*1,NKN0TY 

LEFTY (JJ)*LEFTX(JJ) 

DO 272 KK=1 ,2 
DO 271 I 1*1, NY 

YSPLINEC II , JJ,KK)*XSFLINE( II , JJ,KK) 

271 CONTINUE 

272 CONTINUE 

27 CONTINUE 
GOTO 293 

28 CONTINUE 

DO 29 J=1,NKN0TY 

CALL I NTERV ( T Y , NUHY ,B( J) , LEFTY C J ) , NFL AG > 

CALL BSPLVD ( TY , KY , B < J > , LEFT Y < J > , WORK , DB I ATX , JDERIV) 

DO 292 JJ*1 ,2 

DO 291 11*1, KY 

JB*LEFTY( J)-KYHI 

IF(JB.LE.O) GOTO 291 

TSPLINE<JB,J,JJ)=DBIATX(II,JJ) 

291 CONTINUE 

292 CONTINUE 

29 CONTINUE 

293 CONTINUE 
RETURN 
END 


SUBROUTINE TENVALF ( ARR , LEFTX , LEFTY , KX , KY 
* , I, J, VALUE, I DERI V, JDERIV) 

TENVALF CONFUTES PARTIAL DERIVATIVES FOR A TENSOR 
PRODUCT SPLINE FUNCTION AS INDICATED BY THE PARAHETERS 


00762 


IDERIV, JDERIV. 
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007431C 



00744 :c 

INPUT 


00745 :c 



00744 SC 

ARR<* • 

ARRAY OF COEFFICIENTS FOR SPLINE 

0074? SC 


FUNCTION 

00748 SC 

I » J. * • 

INDICIES FOR POINT OF EVALUATION 

00749.C 


(A(I),B(J)),MCRE A,B ARE ARRAYS 

00770SC 


WHICH CONTAIN COORDINATES OF A 

00771 SC 


SQUARE MESH. 

00772 SC 

LEFTX . » » 

VALUE INDICATING KNOT INTERVAL 

00773SC 


ON WHICH A(I) LIES 

00774 SC 

LEFTY... 

VALUE INDICATING KNOT INTERVAL 

00775SC 


ON WHICH B(J) LIES 

00774 SC 

KX,KY... 

ORDER OF B-SPLINES IN X-DIRECTION, 

00777SC 


Y DIRECTION 

00778 SC 

1DERIV... 

ORDER OF DERIVATIVE DESIRED FOR X 

0077? SC 


DIRECTION 

00780 SC 

JDERIV... 

ORDER OF DERIVATIVE DESIRED FOR Y 

00781 SC 


DIRECTION 

00782 SC 

(IN COMMON ) 


00783 SC 

XSPLINE,YSPLINE. .ARRAYS CONTAINING FUNCTION 

00784SC 


VALUES AND FIRST DERIVATIVES OF B-SPLINES 

00785 SC 


X DIRECTION, Y DIRECTION AT EACH POINT 

00784SC 


GIVEN IN ARRAYS A,B 

00787 SC 



00788 SC 

OUTPUT 


00789SC 



00790SC 

VALUE... 

VALUE OF TENSOR PRODUCT SPLINE 

00791 SC 


OR DERIVATIVE 

00792 SC 



00793 S 

COMMQN/SPL I NES/XSPL I NE < 50 , 1 00 , 2 > , YSPL I NE ( 50 , 1 00 , 2 ) 

00794 S 

DIMENSION ARRUOO'lOO) 

00795 5 

VALUE*0. 


00794 S 

DO 13 JJ*1, 

KY 

00797 S 

JB-LEFTY-KY+JJ 

00798 5 

IF(JB.LE.O) 

GOTO 13 

0079? S 

DO 12 11*1, 

KX 

00600 S 

IB*LEFTX-KX+II 

00801 S 

IF(IB.LE.O) 

GOTO 12 

00802 S 

VALUE*VALUE+ARR< IB,JB)*XSPLINE( IB, I , IDERIVH > 

00803S 

* 9YSPLINE( JB, J, JDERIV-fl ) 

00804 S 

12 CONTINUE 


00805S 

13 CONTINUE 


00804 S 

RETURN 


00807S 

END 


00808 SC 



0080? SC 



00810 S 

SUBROUTINE 

TENSORVAL( ARR,NX , NY , KX , KY , A , B , VALUE 

00811 S 

* ,JDX, JDY) 


00812SC 

THIS ROUTINE CALCULATES THE VALUE OF A TENSOR PRODUCT 

0O813SC 

SPLINE AT THE POINT (A,B)» 


008141C 

00815JC 

oobujc 

00817JC 
00818IC 
00819 JC 

00820 :c 

00821JC 
00822 JC 

00823 :c 

00824!C 
00825 JC 
00824!C 
00827 JC 
00828 1C 
00829 JC 
00830! C 
00831 :c 
00832! 
00833! 
00834! 
00835! 
00834! 
00837! 
00838! 
00839! 
00840! 
00841 !C 
00842 JC 
00843 JC 
00844! 
00845 JC 
00846JC 
00847 !C 
00848 JC 
00849 JC 
00850 JC 
00851 JC 
00852 !C 
00853 JC 
00854 !C 
00855 JC 
00854 !C 
00857 !C 
00858 JC 
00859 JC 
00840 JC 
00841 JC 
00842 JC 
00843 JC 
00844 JC 


INPUT 

ARR. . . ARRAY OF COEFFICIENTS 

NX, NY... DIMENSION OF SPLINE SPACE IN X DIRECT., 

Y DIRECTION. ARRAY ARR WILL HAVE 
DIMENSION NX8NY* 

KX,KY... ORDER OF B-SPLINES IN X DIRECT, Y DIRECT 
A,B... POINT OF EVALUATION 

(IN COMMON ) 

TX,TY. . .KNOT SEQUENCE FOR X DIRECT, 

Y DIRECT. 


OUTPUT 


VALUE... VALUE OF TENSOR PRODUCT SPLINE AT (A,B) 
COMMON/KNOTS/TX ( 100 ) , TY ( 1 00 ) 

DIMENSION BCOEF( 100) ,ARR( 100,100) 

CALL INTERV (TY, NY, B, LEFTY, MFLAG) 

VALUE =0. 

DO 10 J*1 ,KY 

10 BCOEF ( J )*BVALUE ( TX , ARR ( 1 , LEFTY-KY+J > , NX , KX , A , JDX ) 
VALUE*8VALUE ( TY ( LEFTY-KY+1 ) , BCOEF , KY , KY , B , JDY ) 
RETURN 
END 


SUBROUTINE JACOB(NX,NY,KX,KY,A,B) 

THIS ROUTINE COMPUTES THE JACOBIAN OF A 
TENSOR PRODUCT B-SPLINE MAPPING. 


INPUT 

NX, NY... DIMENSION OF SPLINE IN X 
DIRECTION, Y DIRECTION 

KX,KY... ORDER OF B-SPLINES IN X DIRECTION 

Y DIRECTION 

A,B... ARRAYS CONTAINING EVALUATION 

POINTS 

(IN COMMON) 

ALPHA, BETA... COEFFICIENTS OF B-SPLINES 
NKNOTX,NKNOTY. NUMBER OF ELEMENTS IN 
ARRAYS A,B 

OUTPUT < TO TERMINAL) 

A,B," ARRAYS CONTAINING EVALUATION 


00863IC 
00866 :c 
00867 :c 
00868) 
00869) 
00870) 
00871) 
00872) 
00873) 
00874) 
00873) 
00876) 
00877) 
00878) 
00879) 
00880) 
00881) 
00882) 
00883) 
00884) 
00885) 
00886) 
00887) 
00888) 
00889) 
00890) 
00891 )C 
00892 )C 
008931C 
00894) 
00895.C 
00896)C 
00897 )C 
00898 1C 
00899 )C 
00900)0 
00901 )C 
00902)0 
00903)0 
00904)0 
00905)0 
00906)0 
00907)0 
00908 1C 
00909)0 
00910)0 
00911)0 
00912)0 
00913)0 
00914)0 
00915)0 


POINTS 

AJCQBXAN... JACOBIAN AT EACH POINT 

CQHNON/CQEF / ALPHA < 100, 100) ,BETA< 100, 100) 
COMMON/KNOT/NKNOTX , NKNOTY 
REAL A(100) ,B(100) 

PRINT*, ' X Y JACOBIAN' 

PRINT*,' ' 

DO 20 I-1,NKN0TX 
DO 10 J*l, NKNOTY 
II*I 
JJ-J 

CALL TENSORVAL(ALPHA,NX,NY,KX,KY,A(II),B< JJ), 

* XDFIRST, 1,0) 

CALL TENSQRVAL(ALPHA,NX,NY,KX,KY,A(II) ,B(J t), 

* YDFIRST»0,1) 

CALL TENSORVAL < BET A , NX , NY , KX , KY , A < 1 1 ) , B < J J ) , XDSEC , 

* 1 , 0 ) 

CALL TENSORVAL < BET A , NX , NY , KX , KY . A < 1 1 ) , B < J J > , YDSEC , 

* 0 , 1 ) 

AJC0BIAN*XDFIRST*YDSEC-XBSEC*YDFIRST 
PRINT*,A(II),B(JJ),AJCOBIAN 
10 CONTINUE 
20 CONTINUE 
RETURN 
END 


SUBROUTINE CORANGE <NKNOTX f NKNOTY> 

CORANGE DETERMINES THE RANGES OF SUMMATION 
NEEDED TO MINIMIZE THE SMOOTHING FUNCTIONAL G IN 
EACH COORDINATE DIRECTION. 

SPECIFICALLY, FOR EACH COEFFICIENT ALPHA(1,J), 

OR BETA(I, J) , IT DETERMINES THE RANGE OF INDICES FOR 
THE MESH POINTS LYING ON THE SUPPORT OF THE TENSOR 
PRODUCT B-SPLINE HAVING SUBSCRIPTS I,J. IT PLACES 
THE SMALLEST IN IFIRST(I) AND JFIRST(J) AND THE 
LARGEST INDICES IN ILAST(I) AND JLAST(J). THESE VALUES 
DETERMINE WHICH TERMS IN G SHOULD BE SUMMED WHEN 
MINIMIZING THE SMOOTHING FUNCTION IN THE DIRECTION 
REPRESENTED BY THE COEFFICIENT ALPHA(I,J) OR BETA(I,J). 

INPUT 

NKNOTX, NKNOTY... DIMENSIONS FOR SQUARE MESH 

OR NUMBER OF ELEMENTS IN ARRAYS A,E> 

(IN COMMON) 

A,B... ARRAYS CONTAINING COORDINATES FOR 

SQUARE MESH 


00916JC 
00917JC 
00918JC 
00919JC 
00920 JC 

00921 :c 

00922 JC 
00923JC 
00924 JC 
00923 JC 
00926 JC 
00927JC 
00928JC 
00929 JC 
00930 JC 
00931JC 
00932 JC 
00933 JC 
00934 JC 
00935JC 


NX, NY... 


KX,KY... 


LEFTX, LEFTY... 


DIMENSION OF SPLINE SPACE IN X DIRECT 
»Y DIRECTION OR 

TOTAL NO. OF B ■‘SPLINES IN X DIRECTION 

Y DIRECTION 

ORDER OF B-SPLINES IN X DIRECTION, 

Y DIRECTION 

ARRAYS IDENTIFYING KNOT INTERVALS 
ON HHICH SQUARE MESH COORDINATES 
LIE. (LEFTXd)=J IMPLIES 
TX<JX*AdKTX( J+l) 


OUTPUT 

IFIRST , JFIRST . . ARRAYS CONTAINING STARTING 

POINTS FOR THE RANGES OF SUMMATION 
CORRESPONDING TO EACH COEFFICIENT 
ALPHA(I, J) ,BETA(I , J) » 

ILAST, JLAST . » ARRAYS CONTAINING FINAL 

POINTS FOR THE RANGES OF SUMMATION 
CORRESPONDING TO EACH COEFFICIENT 


00936 JC ALPHA(I,J) ,8ETA(I,J> 

00937 JC 

00938 J COMMON/P ARAM2/ A ( 100 ) , B ( 100 ) , NX , NY , KX , K Y , 

00939 J * LEFTX (100), LEFTY (100) 

00940 J C0MM0N/RANGE/IFIRST(100) , ILAST < 100) , JFIRST( 100) , 

00941 J « JLAST (100) 


00942 J 
00943 J 
00944 J 
00945 J 
00946 J 
00947 J 


DO 50 1=1, NX 
IFd.EQ.l) THEN 
11*1 


ELSE 

II*IFIRSTd-l) 

ENDIF 


00948. 10 IFd.GE.LEFTXdl)-(KX-l) .AND, I . LE. LEFTX ( II ) ) 

00949 J * GOTO 20 

00950 J 11*11+1 

00951 J GOTO 10 

00952 J 20 IFIRSTd)*II 

00953 J 30 11*11+1 

00954 J IFdl.GT.NNNOTX) GOTO 40 

00955. IFd.LT.LEFTXdl)-(KX-l)) G n TO 40 

00956 J GOTO 30 

00937 5 40 ILAST < I ) *11—1 

00958 J 50 C1NTINUE 

00939 J DO 100 J=1 ,NY 

00960 J IF(J.EQ.l) THEN 

00961 J JJ*1 

00962 J ELSE 

00963 J JJ*JFIRST< J-l ) 

00964 J ENDIF 

00965 J 60 IF<J,GE,LEFTYUJ)-<KY-1).AND.J.LE.LEFTY<JJ>) 


00966 J * GOTO 70 
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00947? 
00948: 
00949: 
00970: 
00971: 
00972: 
009731 
00974: 
00973: 
00974: 
00977 J 
00978 :c 
00979 :c 
00980 !C 
0098U 
00982 :c 
00983 :c 
00984 :c 
00985 :c 
00984 :c 
00987 :c 
00988 :c 
00989 :c 
00990 :c 
00991 :c 
00992 :c 
00993 :c 
00994 :c 
00995 :c 
00996 :c 
00997 :c 
00998 :c 
00999 :c 
01000JC 

oioonc 
01002 :c 

010031 C 
01004IC 

01005 :c 

01004 JC 

01007: c 
oioo8:c 

oioo9:c 

010105C 


JJ-JJ+1 
GOTO 60 

70 JFIRST<J)-JJ 
80 JJ»JJ+1 

IF(JJ.GT.NKNOTY) GOTO 90 

IF< J*LT«LEFTY< JJ)-(KY-l) ) GOTO 90 

GOTO 80 

90 JLAST( J)*JJ-1 
100 CONTINUE 
RETURN 
END 


REAL FUNCTION GF<II,JJ) 


FUNCTION OF COHPUTES THE SUH OF THE TERMS IN 
THE SMOOTHING FUNCTIONAL 0 OVER THE RANGES INDICATED 
BY IFIRST(II) » JFIRST< JJ) AND ILAST(II), JLAST(JJ). 

INPUT 

II, JJ... INDICES FOR COEFFICIENT INVOLVED 

IN MINIMIZATION. 

(IN COMMON) 

NKNOTX,NKNOTY... DIMENSIONS FOR SQUARE MESH 

OR NUMBER OF ELEMENTS IN ARRAYS A,B 
A, B. . . ARRAYS CONTAINING COORDINATES FOR 

SQUARE MESH 

NX, NY... DIMENSION OF SPLINE SPACE IN X DIRECT. 

, Y DIRECTION OR 

TOTAL NO, OF B-SPLINES IN X DIRECTION, 

Y DIRECTION 

KX,KY. . . ORDER OF B-SPLINES IN X DIRECTION, 

Y DIRECTION 

LEFTX, LEFTY... ARRAYS IDENTIFYING KNOT INTERVALS 
ON WHICH SQUARE MESH COORDINATES 
LIE. < LEFTX ( I ) * J IMPLIES 
T X < J ) < e A ( I ) <TX ( J+ 1 ) 

IFIRST, JFIRST . . ARRAYS CONTAINING STARTING 

POINTS FOR THE RANGES OF SUMMATION 
CORRESPONDING TO EACH COEFFICIENT 


01011 :c ALPHA ( I , J) ,BETA( I , J) . 

01012IC 1LAST, JLAST , « ARRAYS CONTAINING FINAL 

01013:C POINTS FOR THE RANGES OF SUMMATION 


0101V.C CORRESPONDING TO EACH COEFFICIENT 

010151C ALPHA(I,J) ,BETA(I,J) 

01016IC W1,W2... WEIGHTS FOR JACOBIAN, DOT PRODUCT 

oioi7:c 




01018JC 

oioiP! 1 : 

OUTPUT 

O1020JC 

OF... PARTIAL SUM OF TERMS IN G OVER THE 

01021 JC 
01022 !C 


APPROPRIATE RANGE. 

01023! 


REAL AJ<30,30),D0T<30,30> 

01024! 


C0MN0N/PARAM2/A < 100 ) , B ( 1 00 ) , NX , NY , KX , KY , 

01023! 

* 

LEFTXC 100) , LEFTY <100) 

01026! 


COMMON/COEF/ ALPHA < 1 00 , 1 00 ) , BETA < 1 00 , 1 00 ) 

01027! 


COMMON/KNOT /NKNOTX , NKNOTY 

01028! 


C0MM0N/UEIGHTS/U1 ,U2 

01029! 


COMMON/RANGE/IFIRST < 100) , ILAST < 100) , JFIRST dOO) » 

01030! 

* 

JL AST (100) 

01031 ! 


NUMX-NKNOTX 

01032! 


NUMY -NKNOTY 

01033! 


DELX=1 ./(NUMX-1 , ) 

01034! 


DELY*1 ./(NUMY-l . ) 

01035! 


SDELX»DELX*DELX 

01036! 


SDELY*DELY*OELY 

01037! 


SUM-0.0 

01038! 


IF=IFIRST ( II ) 

01039! 


JF-JFIRST ( JJ) 

01040! 


IL-ILAST ( II ) 

01041! 


JLsJLAST(JJ) 

01042! 


IF(IF.GT.l) IF* IF— 1 

01043! 


IF(JF.GT.l) JFaJF-1 

01044! 


IFUL.LT.NUMX) IL-IL+1 

01045! 


IF(JL.LT.NUMY) JL-JL+1 

01046! 


DO 200 J«JF,JL 

01047! 


DO 100 I = IF , IL 

01048! 


CALL TENVALF < AL. HA , LEFTX < I ) , LEFT Y ( J > , KX , KY , I , J , 

01049! 

* 

FIX, 1,0) 

01050! 


CALL TENVALF < BETA, LEFTX < I), LEFTY <J),KX,KY, I, J, 

01051! 

* 

F1Y,1,0) 

01052! 


CALL TENVALF < ALPHA , LEFTX <1 ) , LEFTY ( J ) , KX , N Y , I , J , 

01053! 

* 

F2X,0,1> 

01054! 


CALL TENVALF (BETA. LEFTX ( I > , LEFTY ( J) ,NX»KY, I,J, 

01055! 

* 

F2Y,0, 1 ) 

01056! 


AJ<I,J)=F1X*F2Y-F2X*F1Y 

01057! 


D0T<I,J)=F1X*F2X+F1Y*F2Y 

01058! 


SUM*SUM+D0Td,J)«2 

01059! 

100 

CONTINUE 

01060! 

200 

CONTINUE 

01061! 


SUM1-0. 

01062! 


SUM2-0 . 

01063! 


DO 400 J=JF , JL-1 

01064! 


DO 300 I=IF,IL-1 

01065! 


IF (I .EQ. 1 . AND. J.EQ. 1 ) GOTO 300 

01066! 


IFd.EQ.l. AND. J.EQ. NUMY-l) GOTO 300 

01067! 


IF(I.EQ. NUMX-1. AND. J.EQ.l) GOTO 300 

01068! 


SUM1 *SUM 1 + < A J < I , J ) -A J ( I , J-H ) ) **2 
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01069: 

01070! 

01071: 

01072! 
01073! 
01074! 
01073! 
01076! 
O1077!C 
01078! C 
0107?!C 
01080! 
01081!C 
010821C 
01083JC 
010B41C 
01085! C 
0108AJC 
01087JC 
01088JC 
01089JC 
01090 !C 
oio 9 i:c 
01092JC 
010931C 
010941C 
01093 !C 
0109AJC 
01097 :c 
01098JC 
01099JC 
01100JC 
011011C 
01102JC 
01103JC 
011041C 

oi 105 : c 

01106!C 

01107’C 

0U08JC 

01109JC 

oi i io : c 
oim :c 

0U121C 

01113SC 

011141C 

oiusic 

01U6JC 
011171C 
01U8!C 
Oil 19? C 


SUM2*SUH24 < 6 J ( 1 + 1 , J ) -6 J < I , J ) ) **2 
300 CONTINUE 
400 CONTINUE 

SUH1-SUM14DELX/DELY 

SUH2*SUM2*DEL Y / DEL X 

GF*W1 * < SUN 1 +SUM2 ) + W24DELX *DEL Y *$UH 

RETURN 

END 


SUBROUTINE FFHIN(ERMAX) 


FFNIN SEARCHES FOR THE HIN1HUH OF THE SMOOTHING 
FUNCTIONAL G. EACH CALL TO FFNIN PRODUCES ONE 
CONPLETE ITERATION OF THE CYCLIC COORDINATE METHOD, 
A HULTIDIHENS10NAL SEARCH TECHNIQUE FOR NININIZING 
A FUNCTION OF SEVERAL VARIABLES WITHOUT USING 
DERIVATIVES. THE ROUTINE SEARCHES FOR A MINIMUM 
ALONG EACH COORDINATE DIRECTION. 

IN FFNIN THE COORDINATE DIRECTIONS ARE REPRE- 
SENTED BY THE TENSOR PRODUCT COEFFICIENTS .FOR EACH 
COEFFICIENT THE ROUTINE FIRST DETERHINES THE 
INTERVAL ON WHICH THE COEFFICIENT NUST LIE IF THE 
JACOBIAN OF THE TENSOR PRODUCT HAPPING IS TO BE 
NONNEGATIVE AT ALL NESH POINTS AFFECTED BY THE 
COEFFICIENT. IT THEN CALLS EITHER TESTMINO, TEST- 
MINL, TESTMINR, OR TESTNINB DEPENDING ON WHETHER 
THE INTERVAL IS BIINFINITE, HAS A LEFT ENDPOINT, 

A RIGHT ENDPOINT, OR TWO ENDPOINTS. THE CHOSEN 
SUBROUTINE FINDS THE LOCATION OF THE MININUN OF 
GF ON THE INTERVAL AND CHANGES THE APPROPRIATE 
COEFFICIENT ACCORDINGLY. 


INPUT 


(IN COMMON) 
ALPHA, BET A.. 

A , B » , . 

NX, NY... 

KX,KY » . . 

LEFTX, LEFTY... 


ARRAYS CONTAINING COEFFICIENTS OF 
TENSOR PRODUCT SPLINE MAPPING. 

ARRAYS CONTAINING COORDINATES FOR 
SQUARE NESH. 

DIMENSION OF SPLINE SPACE IN X DIRECT. 
,Y DIRECTION OR 

TOTAL HO. OF B-SPLINES IN X DIRECTION, 

Y DIRECTION. 

ORDER OF B-SPLINES IN X DIRECTION, 

Y DIRECTION. 

ARRAYS IDENTIFYING KNOT INTERVALS 
ON WHICH SQUARE NESH COORDINATES 
LIE. (LEFTX(I)aJ IMPLIES 
TX( JX»A( I )<TX( J+l) ) 


01120* C 



oii2i;c 

1FIRST, JFIRST . » ARRAYS CONTAINING STARTING 

011221c 

POINTS FOR THE RANGES 

OF SUMMATION 

01123* C 

CORRESPONDING TO EACH 

COEFFICIENT 

011241C 

ALPHA! I, J> 'BETA! I, J). 


0U251C 

1LA8T, JLAST « « ARRAYS CONTAINING FINAL 

01126:c 

POINTS FDR THE RANGES OF SUMMATION 

01127JC 

CORRESPONDING TO EACH 

COEFFICIENT 

0U281C 

ALPHAd,J),BETAd,J)» 


011291C 

Ml, M2... HEIGHTS FOR JACOBIAN, 

DOT PRODUCT 

oil 30 : c 

TO BE USED IN GF. 


0U3iic 

NKNOTX,NKNOTY... DIMENSIONS FOR SQUARE 

MESH 

011321C 

0U331C 

OR NUMBER OF ELEMENTS 

IN ARRAYS A,l 

01 134 I C 

OUTPUT 


011351C 

011361C 

ERMAX... MAXIMUM CHANGE IN THE 

COEFFICIENTS 

011371C 

AFTER A COMPLETE ITERATION. 

011381C 

(IN COMMON) 


01 139 • C 

ALPHA, BETA... ARRAYS CONTAINING NEW COEFFICIENTS 

011401C 

FOR TENSOR PRODUCT SPLINE MAPPING. 

01141 SC 
011421 

COMMON/COEF /ALPHA ( 100 , 100 ) , BET A ( 100 , 

100) 

011431 

CQHMQN/PARAM2/A< 100) ,B( ICO) ,NX,NY ,KX 

,KY f 

011441 

* LEFTX < 100 ) , LEFTY ( 100 ) 


011431 

COMHON/RANGE/IFIRST(100),ILAST(100), 

JFIRST ( 100) , 

01146! 

* JLAST (100) 


01147! 

COMMON/PARAM/FKOUNT 


011481 

CQMM0N/UEISHTS/U1 , U2 


01149! 

COMMON/KNOT/ NKNOTX , NKNOTY 


011501 

REAL LEND, LINT, AJ(2) 


011311 

INTEGER CTEST 


011521 

FK0UNT*0 


011531 

ERMAX*0. 


011341 

DO 500 I*2,NX-1 


011531 

DO 400 J*2,NY-1 


011561 

IF*IFIRST ( I ) 


011571 

IL*ILAST(I) 


011581 

JF=JFIRST(J) 


01159! 

JL» JLAST (J) 


011601 

MK0UNT*0 


011611 

5 CONTINUE 


011621 

HOLD*GF < I , J) 


011631 

MTESr=MK0UNT/2*2 


011641 

IF (MTEST .EQ.HKOUNT) THEN 


011651 

HOCO*ALPHA( I , J) 


01166! 

ELSE 


01167! 

HOCO»BETA( I, J) 


01168! 

ENDIF 


01169! 

LEND*-10.0E8 


01170! 

REND-10.0E8 



01171! 

01172! 

01173! 

01174! 

01173! 

01176! 

01177! 

01178! 

0117?! 

01180! 

01181! 

01182! 

ones: 

01184! 

01185! 

01186! 

01187! 

01188! 

01189! 

01190! 

01191! 

01192! 

01193! 

01194! 

01195! 

01196! 

01197! 

01198! 

01199! 

01200 ! 

01201 ! 

01202 ! 

01203! 

01204! 

01203! 

01206! 

01207! 

01208! 

01209! 

01210 ! 

01211 ! 

01212 ! 

01213! 

01214! 

01215! 

01216! 

01217! 

01218! 

01219! 

01220 ! 

01221 ! 


KFLA8*0 
I FI 40*0 

DO 200 II*IF,IL 

DO 100 JJ-JF f JL 

IF ( RKOUNT /2*2 . EQ . HKQUNT > THEN 

ALPHA<I,J)*0, 

ELSE 

BETA I,J>*0. 

ENDIF 

DO IT K-1,2 

CALL TENVALF < ALPHA , LEFTX ( I I ) , LEFT Y ! J J ) , KX , K Y , I I , 

* JJ'FIX'1,0) 

CALL TENVALF < BETA , LEFTX ! 1 1 > , LEFT Y < J J > , KX , KY , 1 1 , 

* JJ,F1Y,1,0) 

CALL TENVALF < ALPHA , LEFTX ( II ) , LEFTY ( J J ) , KX , KY , 1 1 , 

* JJ,F2X,0»1) 

CALL TENVALF! BETA,LEFTX( II ) »LEFTY(JJ) ,KX,KY,II, 
f JJ,F2Y,0,1> 

AJ!K)»F1X*F2Y-F2X*F1Y 
IF!HK0UNT/292*EQ*HK0UNT) THEN 
ALPHA! I , J)*l • 

ELSE 

BETA(I, J)»l. 

ENDIF 

10 CONTINUE 
D*AJ( 1 ) 

C*AJ!2)-D 

IF!C.BT,1.0E-7) THEN 
LINT— D/C 

IF ! LINT . GT. LEND) THEN 
IF < L INT «LEiREND> THEN 
LEND*LINT 
ELSE 
I FLAG*- 1 
ENDIF 
ENDIF 

ELSE IF!C.LT.-1.0E-7> THEN 
RINT—D/C 

IF!RINT.LT.REND) THEN 
IF!RINT .GE.LEND) THEN 
REND*RINT 
ELSE 

IFLAG— 1 
ENDIF 
ENDIF 
ELSE 

KFLAG*KFLAG+1 

ENDIF 

100 CONTINUE 
200 CONTINUE 

CTEST*! ILAST ( I )-IFIRST < I HI X! JLAST! J)-JFIRST( JHl ) 
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01222! IF<KFLAG*EQ*CrEST*OR*IFLAG»EQ»-l> THEN 

01223* IF <HK0UNT/2*2.EQ* MKOUNT) THEN 

01224* ALPHA < I , J ) * < LEND* REND ) /2 ♦ 

01225! ELSE 

01226! BETA<I,J)=(LENDfREND)/2. 

01227! ENDIF 

01228! ELSE IF<LEND.LT.-1 . 0E7* AND. REND. GT .1 ,0E7) THEN 

0122?! CENTER=0.0 

01230! CALL TE5TMIN0<MK0UNT,I »J) 

01231! ELSE IF(LEND.LT.-10.0E7> THEN 

01232! CENTER=0.0 

01233! CALL TESTMINR<MKOUNT,I,J,LEND,REND) 

01234! ELSE IFtREND.GT. 10.E7) THEN ! 

01235! CENTERsO.O * j 

01236! CALL TESTMINL ( MKOUNT , I , J , LEND , REND ) 

01237! ELSE 

01238! CENTER* < LEND+REND ) /2 . 

01239! CALL TESTMINB (MKOUNT, I ,J, LEND, REND) 

01240! ENDIF 

01241! S2*GF( I, J) 

01242! DIFF=H0LD-S2 

01243! IF(HOLD»LT .S2) THEN 

01244! IF (MTEST.EQ. MKOUNT) THEN 

01245! ALPHA ( I , J ) *HOCO 

01246! ELSE 

01247! BETA ( I , J)=HQCO 

01248! ENDIF j 

01249! S2=H0LD I 

01250! DIFF=0. j 

01251! ENDIF 

01252! IF(ABS(DIFF) ,GT .ERMAX) ERMAX=ABS(DIFF) j 

01253! PRINT*, FUNCTION VALUE IS' ,S2 

01254! PRINT*, 'COUNT IS',FKOUNT j 

01255! KFLAG=0 i 

01256! MKOUNT =MKOUNT +1 [. 

01257! IF ( MKOUNT. NE.MK0UNT/2*2) GOTO 5 

01258! 400 CONTINUE 

01259! 500 CONTINUE 

01260! RETURN 

01261! END 

01262IC 

01263!C 

01264JC 

C 1265 ! SUBROUTINE CRIT ( CENTER , C 1 , C2 , C3 , C4 , C5 , NROOTS , R , MKOUNT , 

01266! X I, J) 

012675C 

01268JC CRIT FINDS THE COEFFICIENTS OF THE 4TH DEGREE 

012695C POLYNOMIAL REPRESENTING GF AND COMPUTES ITS CRITICAL 
01270JC POINTS, I *E. , IT FINDS THE POINTS FOR WHICH THE DERIVATIVE 
01271 !C OF THE POLYNOMIAL IS 0. 

012721C 
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012731 C 
012741C 
01275 1C 
012761C 
012771C 
012781C 
012791C 
01280 1C 
01281 tC 
01282 !C 
012831C 
01284JC 

01285 :c 

012861C 

01287 ;c 
0128810 
01289 1C 
01290 1C 
01291 1C 
01292 1C 
012931C 
012941C 
01295 JC 
01296 1C 
01297JC 
01298 1C 
01299 1C 
01300 I C 
013011C 
01302 1C 
01303 1C 
013041C 
01305 1C 
013061C 
01307 1C 
01 308 1C 
013091C 
01310 1C 
013111C 
0131 2 1 C 
013131C 
013141C 
01315 1 C 
013161C 
013171C 
013181C 
013191 
013201 
013211 
013221 
013231 


INPUT 


CENTER** ♦ 


HKOUNT * . . 


I,J*. 

(IN COMMON) 
ALPHA, BETA. . 

AiBi ♦ ♦ 

NX, NY... 

KXrKY . * * 

LEFTX, LEFTY * * . 


NUMBER AT CENTER OF INTERVAL TO BE 
CONSIDERED. IF INTERVAL IS INFINITE THEN 
CENTER ASSIGNED A VALUE OF 0* 

HKOUNT EVEN MEANS THE COEFFICIENT 
INVOLVED IN MINIMIZATION IS IN THE ALPHA 
ARRAY. MKOUNT ODD MEANS THE COEFFICIENT 
IS IN THE BETA ARRAY. 

SUBSCRIPTS FOR COEFFICIENT INVOLVED 
IN MINIMIZATION 

ARRAYS CONTAINING COEFFICIENTS OF 
TENSOR PRODUCT SPLINE MAPPING. 

ARRAYS CONTAINING COORDINATES FOR 
SQUARE MESH. 

DIMENSION OF SPLINE IN X DIRECTION 
,Y DIRECTION OR 

TOTAL NO. OF B-SPLINES IN X DIRECTION, 

Y DIRECTION. 

ORDER OF B-SPLINES IN X DIRECTION, ' 

Y DIRECTION. 

ARRAYS IDENTIFYING KNOT INTERVALS 
IN WHICH SQUARE MESH COORDINATES 
LIE. (LEFTX(I)=J IMPLIES 
TX( JK=A(IKTX( J+l) ) 


U1,U2... WEIGHTS FOR JACOBIAN, DOT PRODUCT 

TO BE USED IN GF. 

NKNOTX,NKNOTY. ..DIMENSIONS FOR SQUARE MESH 

OR NUMBER OF ELEMENTS IN ARRAYS A,B. 

FKOUNT • * , PARAMETER CONTAINING NUMBER OF CALLS TO 

GF 

XK... ARRAY CONTAININ POINTS -2, -1,0, 1,2 

WHIC“ ARE USED AS TEST POINTS IN 
DETERMINING THE COEFFICIENTS OF THE 
4TH DEGREE POLYNOMIAL WHICH APPROXI- 
MATES GF, 


OUTPUT 

C1,C2,C3,C4,C5» « COEFFICIENTS OF 4TH DEGREE POLYNOMIAL. 

Cl IS THE COEFFICIENT OF THE 4TH DEGREE TERN. 
NROOTS... NUMBER OF CRITICAL POINTS 

R... ARRAY CONTAINING CRITICAL POINTS 

REAL R(3),BK(5) 

COMMON/COEF/ALPHA (100,100), BETA (100,100) 

COMMON/ P ARAM2/A (100),B(100),NX,NY,KX,KY,LEFTX(100) 

* , LEFTY ( 100) 

COMMON/KNOT /NKNOTX » NKNOTY 


01324! 
01325! 
01326! 
01327! 
01328! 
01329! 
01330! 
01331! 
01332! 
01333! 
01334! 
01335! 
01336! 
01337! 
01338! 
01339! 
01340! 
01341! 
01342! 
01343! 
01344! 
01345! 
01346! 
01347! 
.01348! 
01349! 
01350! 
01351! 
01352! 
01353! 
01354! 
01355! 
01356! 
01357! 
01358! 
01359! 
01360! 
01361! 
01362!C 
01363!C 
01364!C 
01365! 
01366!C 
01367!C 
01368! C 
01369JC 
01 370 : c 


C0MM0N/WEI6HTS/U1 ,U2 
CONMQN/PARAM/FKOUNT 
COHHON/XKS/XK ( 5 ) 

I F < MK0UNT/2#2 • EQ 4 MKQUNT ) THEN 

DO 100 IK- 1,5 

ALPHA < I, J)-XK(IK) 4CENTER 

BK(IK)*GF(I,J) 

FKQUNT-FKOUNT+1 
100 CONTINUE 
ELSE 

DO 200 IK=1,5 

BET A ( I , J > *XK < IK ) +CENTER 

BK( IK)=GF< I , J) 

FKOUNT *FKOUNT 4 1 
200 CONTINUE 
ENDIF 
D S XK(4> 

B1=BK( 1 ) 

B2=BK<2> 

B3=BK(3) 

B4*BK(4) 

B5=BK(5> 

C5=B3 

SUM=-B5+8.*(B4-B2)+B1 

C4=1./(12.*D)*SUM 

SUM=-B5+16.t(B4+B2>-30.*B3-Bl 

C3=1./(24.*D*D)*SUM 

SUH*B5-2,t(B4-P2)-Bl 

C2=1./(12.*D*D*D)*SUM 

SUM=B5-4.*(B4+B2)+6.*B3+B1 

C1=1./(24.#D**4)*SUM 

IF ( ABS(C1 ) .LT. I .0E-06) GOTO 300 

CALL CUB I C ( 4 . *C1 , 3 . *C2 , 2 . *C3 , C4 , NROOTS , R ) 

RETURN 

300 CONTINUE 
NR00TS=-1 
RETURN 
END 


SUBROUTINE T£STMINO(MNOUNT , I , J) 

FOR A GIVEN COEFFICIENT ALPHA(I,J) OR BETA(I,J> 
TESTMINO FINDS AND TESTS THE CRITICAL POINTS OF 
THE 4TH DEGREE POLYNOMIAL REPRESENTING GF TO 
DETERMINE WHICH POINT YIELDS THE SMALLEST VALUE FOR GF 


01371 !C WHEN GF IS VIEWED AS A FUNCTION OF THAT COEFFICIENT. 
01372IC THE NUMBER CHOSEN BECOMES THE NEW VALUE FOR 
01373IC ALPHA ( I, J) OR BETA(I,J) . 

01374!C 
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013735C 

INPUT 


01 376* C 



013771C 

MKOUNT . * . 

MKOUNT EVEN MEANS THE COEFFICIENT 

01378JC 


INVOLVED IN MINIMIZATION IS IN THE ALPHA 

01379JC 


ARRAY. MKOUNT ODD MEANS THE COEFFICIENT 

01 380 5 C 


IS IN THE BETA ARRAY. 

01381 :c 

I,J.. 

SUBSCRIPTS FOR COEFFICIENT INVOLVED 

01382JC 


IN MINIMIZATION 

01383JC 

(IN COMMON) 


01384 :c 

ALPHA, BETA*. 

ARRAYS CONTAINING COEFFICIENTS OF 

01385JC 


TENSOR PRODUCT SPLINE MAPPING. 

01386 tC 

A , B » » » 

ARRAYS CONTAINING COORDINATES FOR 

01387 »C 


SQUARE MESH. 

01388JC 

NX, NY... 

DIMENSION OF SPLINE IN X DIRECTION 

01389JC 


,Y DIRECTION OR 

01390 ;c 


TOTAL NO. OF B-SPLINES IN X DIRECTION, 

01391 :c 


Y DIRECTION. 

01392 :c 

KX,KY. . . 

ORDER OF B-SPLINES IN X DIRECTION, 

01393IC 


Y DIRECTION. 

01394JC 

LEFTX, LEFTY... 

ARRAYS IDENTIFYING KNOT INTERVALS 

01395J C 


IN WHICH SQUARE MESH COORDINATES 

01396JC 


LIE. (LEFTX( I )=J IMPLIES 

01397JC 


TX< J)<=A(I KTX( J+l) ) 

01398JC 



01399JC 

U1,U2... 

WEIGHTS FOR JACOBIAN, DOT PRODUCT 

014001C 


TO BE USED IN GF. 

01401 :c 

NKNQTX.NKNOTY . 

..DIMENSIONS FOR SQUARE MESH 

01402:c 


OR NUMBER OF ELEMENTS IN ARRAYS A,B. 

01403JC 

FKOUNT . * . 

PARAMETER CONTAINING NUMBER OF CALLS TO 

01404JC 


GF 

01403JC 

XK... 

ARRAY CONTAINING POINTS -2, -1,0, 1,2 

01 406 1C 


WHICH ARE USED AS TEST POINTS IN 

01407JC 


DETERMINING THE COEFFICIENTS OF THE 

01408 JC 


4TH DEGREE POLYNOMIAL WHICH APPROXI- 

01409JC 


MATES GF. 

01410JC 



0141 1 : c 

OUTPUT 


01412SC 



01 4 1 3 : c 

ALPHA(I, J) OR 

BETA( I , J) . .NEW VALUE FOR COEFFICIENT 

01 4 1 4 : c 



01415: 

COMMON/COEF/ ALPHA < 100 , 1 00 ) , BETA ( 100 , 1 00 ) 

014165 

C0MM0N/PARAM2/A ( 1 00 ) , B ( 100 ) , NX , NY , KX , KY , LEFTX (100), 

01417: 

* LEFTY (100) 


01418: 

COMMON/KNOT/NKNOTX ,NKNOTY 

01419: 

C0MM0N/WEIGHTS/U1 ,U2 

01420: 

COMHON/PARAM/FKOUNT 

014211 

COMMON/XKS/XK ( 5 ) 

01422; 

REAL R(3) 


01423: 

FM ( R ) =C1 *R**4+C2*R*R*R+C3*R*RK4»R+C5 

01424: 

BO 30 IK=1,5 


01425: 

XK ( IK ) =FLOAT ( IK) -3 . 
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01426: 

01427J 

01428: 

01429: 

01430: 

01431: 

01432! 

01433! 

01434! 

01435! 

01436: 

01437: 

01438: 

01439: 

01440: 

01441! 

01442! 

01443: 

01444! 

01445! 

01446JC 

01447!C 

01448’C 

01449! 

01450!C 

01451!C 

01452!C 

01453!C 

014541C 

01455!C 

01456JC 

01457!C 

014581C 

01459!C 

01460!C 

01461!C 

01462!C 

01463:c 

01464!C 

01465JC 

01466JC 


50 CONTINUE 

CALL CRIT ( 0 , Cl , C2 , C3, C4 , C5 , NROOTS , R , MKOUNT , I , J ) 

IF < NROOTS. NE.-l) GOTO 55 

RETURN 

55 CONTINUE 
TNIN=10.0E10 
DO 600 IR*1, NROOTS 
FMINN=FM(R(IR)) 

IF(FMINN.LT.TMIN) THEN 

TNIN=FMINN 

IMIN*IR 

ENDIF 

600 CONTINUE 

IF ( MK0UNT/2*2 . EQ . MKOUNT ) THEN 
ALPHA ( I , J) S R(IMIN) 

ELSE 

BETA(I » J)=R(IMIN) 

ENDIF 

RETURN 

END 


SUBROUTINE TESTMINR< MKOUNT . I , J, LEND , REND > 

FOR A GIVEN COEFFICIENT AIPHA(I,J) OR BETA(I,J) 
TESTMINR FINDS AND TESTS THE CRITICAL POINTS OF 
THE 4TH DEGREE POLYNOMIAL REPRESENTING GF TO 
DETERMINE WHICH POINT YIELDS THE SMALLEST VALUE FOR GF 
WHEN GF IS VIEWED AS A FUNCTION OF THAT COEFFICIENT. 

THE SMALLEST VALUE IS COMPARED WITH THE VALUE AT THE 
RIGHT ENDPOINT OF THE INTERVAL (LEND, REND) 

TO DETERMINE AT WHAT NUMBER THE MINIMUM VALUE 
OF GF OCCURS. THE NUMBER CHOSEN BECOMES THE NEW 
VALUE FOR ALPHA ( I, J) OR BETA(I,J) . 

INPUT 

MKOUNT... MKOUNT EVEN MEANS THE COEFFICIENT 

INVOLVED IN MINIMIZATION IS IN THE ALPHA 
ARRAY, MKOUNT ODD MEANS THE COEFFICIENT 


01467!C 
01468JC 
01469IC 
014701C 
01471 :c 
01472!C 


If J.. 


LEND, REND.. 


IS IN THE BETA ARRAY. 

SUBSCRIPTS FOR COEFFICIENT INVOLVED 

IN MINIMIZATION 

LEFT AND RIGHT ENDPOINTS FOR 

INTERVAL. LEND IS A NEGATIVE NO. WITH VERY 

LARGE MAGNITUDE INDICATING THAT THE LEFT 


01473IC 

01474IC 

01475.C 

01476!C 


ENDPOINT IS INFINITE. 

(IN COMMON) 

ALPHA, BETA.. ARRAYS CONTAINING COEFFICIENTS OF 
TENSOR PRODUCT SPLINE MAPPING, 


01477*.C 

ArB* * * 

ARRAYS CONTAINING COORDINATES FOR 

01478JC 


SQUARE MESH. 

01479JC 

NX, NY... 

DIMENSION OF SPLINE IN X DIRECTION 

01480.C 


,Y DIRECTION OR 

01481 iC 


TOTAL NO. OF B-SPLINES IN X DIRECTION 

01482*0 


Y DIRECTION. 

01483? C 

KX,KY.. • 

ORDER OF 8-SPLINES IN X DIRECTION, 

01484*0 


Y DIRECTION. 

01485.*C 

LEFTX, LEFTY... 

ARRAYS IDENTIFYING KNOT INTERVALS 

01486:0 


IN WHICH SQUARE MESH COORDINATES 

01487JC 


LIE. <LEFTX(I)=J IMPLIES 

01488:0 


TX(J)<*A<IXTX< J+l) ) 

01489JC 



01490:0 

W1 r U2. . . 

WEIGHTS FOR JACOBIAN, DOT PRODUCT 

0149i:c 


TO BE USED IN GF. 

01492*0 

NKNOTX , NKNOTY . 

..DIMENSIONS FOR SQUARE MESH 

01493*0 


OR NUMBER OF ELEMENTS IN ARRAYS A, 8. 

01494*0 

FKOUNT * ♦ » 

PARAMETER CONTAINING NUMBER OF CALLS 

01495JC 


GF 

01496:0 

XK ... 

ARRAY CONTAININ POINTS -2, -1,0, 1,2 

01497:0 


WHICH ARE USED AS TEST POINTS IN 

01498:0 


DETERMINING THE COEFFICIENTS OF THE 

01499:0 


4TH DEGREE POLYNOMIAL WHICH APPROXI- 

oi30o:c 


MATES GF. 

oi30i:c 



01302!C 

OUTPUT 


01503*0 



01504JC 

ALPHA ( I , J) OR 

BETA( I , J) • .NEW VALUE FOR COEFFICIENT 

01505?C 



01506: 

COMMON/COEF/ALPHA( 100, ICO) ,Bt~A< 100, 100) 

01507: 

COMMON/PARAM/FKOUNT 

01508: 

C0MM0N/PARAM2/A ( 100 ) , B < 100) , NX , NY , KX , K Y 

01509: 

* ,LEFTX< 100) , 

LEFTY(IOO) 

oi5io: 

COMMON/KNOT /NKNOTX , NKNOTY 

01511,* 

C0MM0N/WEIGHTS/W1 , W2 

01512: 

COMMON/XKS/XK < 5 ) 

01513: 

REAL R(3),LEND 

01514! 

FM<R)»C1*R**4+C2*R*R*R+C3*R*R+C4*R+C3 

01515: 

XK(1)*REND 


01516: 

DO 50 IK=2,5 

. 

01517: 

XK (IK) “REND- 

FLOAT < IK ) 4-1 . 

01518! 

30 CONTINUE 


0151 0 : 

CALL CRIT<0, 

C1,C2,C3,C4,C3,NROOTS,R,MKOUNT,I,J) 

01520! 

IF(NROOTS.NE 

.-1) GOTO 35 

01521*. 

PRINT*,' WARNING *1 COEF IS 0' 

01522: 

RETURN 


01323! 

53 CONTINUE 


01324! 

TMIN=10.E10 


01325! 

IR=0 

1 

01526.* 

DO 600 I ROOT 

*1,NR00T5 

01527.* 

IF(R(IROOT) .LE.REND) THEN 


01528! 

01529! 

01530! 

01531! 

01532! 

01533! 

01534! 

01535! 

01536! 

01537! 

01538! 

01539! 

01540! 

01541! 

01542! 

01543! 

01544! 

01545! 

01546! 

01547! 

01548! 

01549! 

01550! 

01551! 

01552! 

01553! 

01554! 

01555! 

01556! 

01557! 

01558! 

01559! 

01560! 

01561 !C 

01562!C 

01563JC 

01564! 

01565!C 

01566JC 

01567!C 

01568JC 

01569 :c 

01570!C 

01571 !C 

0t572!C 

01573JC 

01574!C 

01575JC 

015765C 

01577IC 

01578!C 


IR*IRF1 

R(IR)=R*IROOT) 

ENDIF 

600 CONTINUE 
NROOTS*IR 

IF(NROOTS.EQ.O) THEN 

IF(MK0UNT/2*2.EG.MK0UNT) THEN 

ALPHA<I,J)=REND 

TMIN*FM(REND) 

ELSE 

BETA(I,J)*REND 

TMIN=FM(REND> 

ENDIF 

ELSE 

DO 700 I ROOT =1 »NROOTS 
FMINN*FM(R( IROOT) ) 
IF*FMINN.LT.TMIN) THEN 
TMIN=FMINN 
IMIN=IROOT 
ENDIF 

700 CONTINUE 

FMINN=FM*REND) 

IF(FMINN.LT.TMIN) R< IMIN)=RENI' 
IF ( MK0UNT/2*2 . EQ . HKOUNT ) THEN 
ALPHA(I,J)=R*IMIN) 

TMIN=FM(R( IMIN) ) 

ELSE 

BETA* I » J)=R( I MIN) 
TMIN=FM(R(IMIN) ) 

ENDIF 

ENDIF 

RETURN 

END 


SUBROUTINE TESTMINL < MKOUNT , I , J, LEND , REND ) 

FOR A GIVEN COEFFICIENT ALPHA* I,J> OR BETA<I,J) 
TESTHINL FINDS AND TESTS THE CRITICAL POINTS OF 
THE 4TH DEGREE POLYNOMIAL REPRESENTING GF TO 
DETERMINE WHICH POINT YIELDS THE SMALLEST VALUE FOR GF 
WHEN GF IS VIEWED AS A FUNCTION OF THAT COEFFICIENT. 
THE SMALLEST VALUE IS COMPARED WITH THE VALUE AT THE 
LEFT ENDPOINT OF THE INTERVAL (LEND, RENIO 
TO DETERMINE AT WHAT NUMBER THE MINIMUM VALUE 
OF GF OCCURS. THE NUMBER CHOSEN BECOMES THE NEW 
VALUE FOR ALPHA* I, J) OR BETA* I, J) . 

INPUT 
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01579JC 

MKOUNT . . * 

MKOUNT EVEN MEANS THE COEFFICIENT 

01590 : c 


INVOLVED IN MINIMIZATION IS IN THE ALPHA 

01581 :c 


ARRAY. MKOUNT ODD MEANS THE COEFFICIENT 

01582 :c 


IS IN THE BETA ARRAY. 

01583 : c 

Z , J * ♦ 

SUBSCRIPTS FOR COEFFICIENT INVOLVED 

01584JC 


IN MINIMIZATION 

01585 JC 

LEND, REND.. 

LEFT AND RIGHT ENDPOINTS FOR 

015865C 


INTERVAL. REND IS A VERY LARGE NUMBER, 

01587JC 


INDICATING THAT THE RIGHT ENDPOINT IS 

01588JC 


INFINITE. 

01589 ;c 

< IN COMMON) 


01590JC 

ALPHA,BETA. . 

ARRAYS CONTAINING COEFFICIENTS OF 

01591 JC 


TENSOR PRODUCT SPLINE MAPPING. 

01592JC 

A ,B . . . 

ARRAYS CONTAINING COORDINATES FOR 

01593JC 


SOUARE MESH. 

01594JC 

NX,NY. .. 

DIMENSION OF SPLINE IN X DIRECTION 

01595JC 


,Y DIRECTION OR 

01596JC 


TOTAL NO. OF B-SPLINES IN X DIRECTION, 

01597JC 


Y DIRECTION. 

01598JC 

KX,KY. . , 

ORDER OF B-SPLINES IN X DIRECTION, 

01599JC 


Y DIRECTION. 

01600JC 

LEFTX, LEFTY. ♦ . 

ARRAYS IDENTIFYING KNOT INTERVALS 

01601 JC 


IN WHICH SQUARE MESH COORDINATES 

01602JC 


LIE, ! LEFTX! I )=J IMPLIES 

01603 JC 


TX( JK=A(I ><TX( J+l) ) 

01604 JC 



01605JC 

U1,U2. , . 

WEIGHTS FOR JACOBIAN, DOT PRODUCT 

01606JC 


TO BE USED IN GF. 

01607JC 

NKNOTX,NKNOTY. 

,, DIMENSIONS FOR SQUARE MESH 

01608 JC 


OR NUMBER OF ELEMENTS IN ARRAYS A.B. 

01609JC 

FKOUNT . . . 

PARAMETER CONTAINING NUMBER OF CALLS TO 

01610 JC 


GF 

01611JC 

XK... 

ARRAY CONTAININ POINTS -2, -1,0, 1,2 

01612JC 


WHICH ARE USED AS TEST POINTS IN 

01613 JC 


DETERMINING THE COEFFICIENTS OF THE 

01614JC 


4TH DEGREE POLYNOMIAL WHICH APPROXI- 

01615 JC 


MATES GF. 

0161 6 J C 



01617 JC 

OUTPUT 


01618JC 



01619JC 

ALPHA! I, J) OR 

BETA! I, J).. NEW VALU" FOR COEFFICIENT 

01620JC 



01621! 

COMMCN/COEF/ALPHA< 100, 100), BETA! 100,100) 

01622 J 

COMMON/PARAM/FKQUNT 

01623! 

C0MM0N/PARAM2/A! 100 ) , B ! 1 00 ) , NX , N Y , KX , K Y 

01624 J 

* , LEFTX! 100) , 

LEFTY! 100) 

01625J 

COMMON/KNOT /NKNOTX , NKNOTY 

01626J 

C0MM0N/UEIGHTS/W1 ,U2 

016275 

C0MM0N/XKS/XK!5> 

01628! 

REAL R!3),LEND 

01629! 

FM! R) «C1*R**4+C2*R*R*R+C3*R*R+C4*R+C5 


01630: 
01631: 
01632: 
01633: 
01634: 
01635: 
01636: 
01637: 
01638: 
01639: 
01640: 
01641: 
01642: 
01643: 
01644: 
01645: 
01646: 
01647: 
01648: 
01649: 
01650: 
01651: 
01652: 
01653: 
01654: 
01655: 
01656: 
01657: 
01658: 
01659: 
01660: 
01661: 
01662: 
01663: 
01664: 
01665: 
01666: 
01667: 
01668: 
01669*. 
01670: 
01671: 
01672: 
01673: 
01674! 
01675: 
01676! 
01677: 
01678! I 
01679:1 
01680: 


XK( 1 )»LEND 
DO 50 IK*2,5 
XK(IK)*LEND+FL0AT(IK>-1 . 

50 CONTINUE 

CALL CRIT ( 0 , Cl » C2 , C3 , C4 » C5 » NROOTS , R , HKOUNT , 
% I,J> 

IF (NROOTS. NE.-l) GOTO 55 
PRINT*, 'WARNING #1 COEF IS 0' 

RETURN 

55 CONTINUE 
THIN*10.0E10 
IL=0 

DO 600 IR=1, NROOTS 
IF ( R ( IR ) . GE . LEND ) THEN 
IL=IL+1 
R(IL)=R(IR) 

ENDIF 

600 CONTINUE 
NROOTS=IL 

IF (NROOTS. EG. 0) THEN 
IF (HKQUNT/2*2.EQ. HKOUNT) THEN 
ALPHA(I , J)=LEND 
THIN=FM(LEND> 

ELSE 

BETA < I . J)=LEND 
THIN- FH( LEND) 

ENDIF 

ELSE 

DO 700 IR=1 ,NROOTS 
FHINN=FH<R( IR) ) 

IF(FHINN.LT.THIN) THEN 

THIN=FHINN 

IHIN=IR 

ENDIF 

700 CONTINUE 

FHINN=FH(LEND) 

IF(FHINN.LT.THIN) R(IHIN)=LEND 
IF(HK0UNT/2*2.EQ. HKOUNT) THEN 
ALPHA ( I » J)=R< IhIN) 

THIN=FH<R<IHIN) ) 

ELSE 

BETA( I, J)=R( IHIN) 

THIN=FH(R(IHIN) ) 

ENDIF 

ENDIF 

RETURN 

END 


SUBROUTINE TESTMINB ( HKOUNT , I , J , LEND , REND ) 
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01681 1C 
016821C 
01683 JC 
01684JC 
01683 JC 
016B6JC 
01687 JC 
01688JC 
0163? JC 
01690JC 
01691 !C 
01692JC 
01693JC 
01694JC 
01693JC 
01696JC 
01697JC 
01698 :c 
01699JC 
01700JC 
01701 :c 
01702JC 
01703JC 
01704JC 
01703JC 
01706JC 
01707JC 
01708SC 
01709JC 
01710*0 
0171 UC 
01712*0 
0171310 
01714JC 
01715*0 


FOR 6 GIVEN COEFFICIENT ALPHA! I, J) OR BETA!I,J> 
TESTHIN8 FINDS AND TESTS THE CRITICAL POINTS OF 
THE 4TH DEGREE POLYNOHIAL REPRESENTING GF TO 
DETERMINE WHICH POINT YIELDS THE SMALLEST VALUE FOR GF 
WHEN GF IS VIEWED AS A FUNCTION OF THAT COEFFICIENT. 
THE SMALLEST VALUE IS COMPARED WITH THE VALUE 
AT THE ENDPOINTS OF THE INTERVAL ( LEND , REND) 

TO DETERMINE AT WHAT NUMBER THE MINIMUM VALUE 
OF GF OCCURS* THE NUMBER CHOSEN BECOMES THE NEW 
VALUE FOR ALPHA! I, J) OR BETA! I, J) . 

INPUT 


MKOUNT 


♦ * ♦ 


I * J* « 

LEND, REND.* 

!IN COMMON) 
ALPHA, BETA*, 

A , B * * * 

NX, NY... 

KX,KY* * * 

LEFTX, LEFTY... 


MKOUNT EVEN MEANS THE COEFFICIENT 
INVOLVED IN MINIMIZATION IS IN THE ALPHA 
ARRAY. MKOUNT ODD MEANS THE COEFFICIENT 
IS IN THE BETA ARRAY, 

SUBSCRIPTS FOR COEFFICIENT INVOLVED 

IN MINIMIZATION 

LEFT AND RIGHT ENDPOINTS FOR 

INTERVAL 

ARRAYS CONTAINING COEFFICIENTS OF 
TENSOR PRODUCT SPLINE MAPPING* 

ARRAYS CONTAINING COORDINATES FOR 
SQUARE MESH. 

DIMENSION OF SPLINE IN X DIRECTION 
»Y DIRECTION OR 

TOTAL NO. OF B-SPLINES IN X DIRECTION, 

Y DIRECTION. 

ORDER OF B-SPLINES IN X DIRECTION, 

Y DIRECTION. 

ARRAYS IDENTIFYING KNOT INTERVALS 
IN WHICH SQUARE MESH COORDINATES 


01716JC 
01717JC 
01718 : C 
01719JC 
01720 JC 
01721JC 
01722JC 
01723 JC 
01724JC 
01725JC 
01726.C 
01727 JC 
01728.C 
01729JC 
01730JC 
01731 JC 


LIE. (LEFTX! I )»J IMPLIES 
TX! JX*A!IXTX( J+l ) ) 

W1,W2* • « WEIGHTS FOR JACOBIAN, DOT PRODUCT 

TO BE USED IN GF. 

NKNOTX,NKNOTY... DIMENSIONS FOR SQUARE MESH 

OR NUMBER OF ELEMENTS IN ARRAYS A,B. 

FKOUNT ♦ • , PARAMETER CONTAINING NUMBER OF CALLS TO 

GF 

XK... ARRAY CONTAININ POINTS -2, -1,0, 1,2 

WHICH ARE USED AS TEST POINTS IN 
DETERMINING THE COEFFICIENTS OF THE 
4TH DEGREE POLYNOMIAL WHICH APPROXI- 
MATES GF. 

OUTPUT 
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01732JC 

01733SC 

01734JC 

01735J 

01736: 

01737J 

01738J 

0173?: 

01740: 

01741 : 

01742: 

01743: 

01744J 

01745: 

01746: 

01747: 

01748: 

01749: 

01750: 

01751*. 

01752: 

01753: 

01754: 

01755: 

01756: 

01757*. 

01758! 

01759*. 

01760: 

01761: 

01762: 

01763: 

01764: 

01765: 

01766: 

01767*. 

01768: 

0176?: 

01770: 

01771: 

01772: 

01773*. 

01774.* 

01775: 

01776: 

01777: 

01778*. 

0177?: 

017805 

01781*. 

01782: 


ALPHAd,J) OR BETAd, J) . .NEW VALUE FOR COEFFICIENT 

COMMON/COEF/ ALPHA ( 100,100) , BETA ( 100, 100) 
COHHON/PARAM/FKQUNT 

C0MM0N/PARAM2/ A ( 1 00 ) , B ( 100 ) , NX , NY , KX , K Y 

* , LEFTX( 100), LEFTY dOO) 

COMMON/KNOT /NKNOTX , NKNOTY 
C0HMQN/WEIGHTS/W1 ,U2 
COMMON/XKS/XK ( 5 ) 

REAL R<3)»LEND 

FH ( R ) *C 1 *R**4+C2*R*R*R+C3*R*R4-C4*R+C5 
CENTER* ( LEND4-RENB ) /2 . 

XKd)=LENB-CENTER 
XK ( 2 ) * ( LEND-CENTER ) /2 . 

XK(3)*0. 

XK < 4 ) * ( REND-CENTER ) /2 . 

XK ( 5 ) *REND-CENTER 

CALL CRI T ( CENTER , C 1 , C2 , C3 , C4 , C5 , NROOTS , R , MKOUNT , 

* I,J) 

IF ( NROOTS. NE.-l) GOTO 55 
RETURN 

55 CONTINUE 

TMIN=10.0E10 

IB*0 

BO 600 I R= UNROOTS 

IF ( RdR ) 4-CENTER . GE . LEND . AND . R < IR ) 4-CENTER . LE . REND ) THEN 

IB*IB+1 

R(IB)=R(IR) 

ENBIF 

600 CONTINUE 
NROOTS* IB 

IF ( NROOTS. EQ.O) THEN 

IF (MK0UNT/242.EQ.MK0UNT) THEN 
ALPHA (I , J ) =CENTER 
THIN*FM(CENTER) 

ELSE 

BETAd, J)*CENTER 
TMIN*FM< CENTER) 

ENBIF 

ELSE 

BO 700 IR* UNROOTS 
FMINN=FM(R( IR) ) 

IF <FMINN.LT. THIN) THEN 

TMIN=FHINN 

IHIN*IR 

ENBIF 

700 CONTINUE 

FMINNL=FH(L£NB) 

FMINNR*FM(RENB> 

IF(FHINNL.LT.FMINNR) THEN 


01783$ 
01784$ 
01785$ 
01786$ 
01787$ 
01788$ 
01789$ 
01790$ 
01791$ 
01792$ 
01793$ 
01794$ 
01795$ 
01796$ 
01797$ 
01798$ 
01799$C 
01800 $C 
01801$C 
01802$ 
01803$C 
01804$C 
01805JC 
01806JC 
01807*0 
01808 $C 
018095C 
01810JC 
01811 JC 
018 12 $ C 
01813$C 
01814JC 
01815SC 
01816JC 
01817JC 
01B18JC 
01819$ 
01820$ 
01821$ 
01822$ 
01823$ 
01824$ 
01825$ 
01826$ 
01827$ 
01828$ 
01829$ 
01830$ 
01831$ 
01832$ 
01833$ 


1F<FHINNL.LT.TMIN> THEN 

R<IMIN)*LEND 

TNIN-FM(LEND) 

EN01F 

ELSE lF(FMINNR.LT.TMlN) THEN 
RUHIN)«REND 
THIN-FH(REND) 

ENDIF 

IF < HKOUNT /2*2 . EQ . MKOUNT ) THEN 
ALPHA< I , J)*RUMIN) +CENTER 
ELSE 

BET A < I , J ) *R < IM IN ) +CENTER 
ENDIF 
ENDIF 
RETURN 
END 


SUBROUT 1 NE CUB I C < A3 , A2 , A 1 , AO , NROOT S , RR > 

CUBIC CONFUTES THE ROOTS OF ft CUBIC POLY- 
NOMIAL USING FORMULAS FROM ‘HANDBOOK OF 
MATHEMATICAL TABLES AND FORMULAS* BY RICHARD 
STEVENS BURINGT0N,PH,D., HCGRAW-HILL ,NEU YORK, 
1962. 


INPUT 

A3 , A2 , A1 , AO . » 


OUTPUT 

NROOTS • » . 

RR ... 

REAL RR(3) 

PI*3.1415 
P-A2/A3 
Q-A1/A3 
R»A0/A3 
A»l./3.*(3.*Q-PtP) 

B»1 . /27 . * < 2 . *P*PtP-9 . tP*Q+27 . *R ) 
IF ( ABS ( B ) . LT . 1 . 0E-06 ) THEN 
SIGNB*0* 

ELSE 

SIGNB>B/ABS(B) 

ENDIF 

BB»B*B/4. 

AAA«A*A*A/27. 

TEST-BB+AAA 


COEFFICIENTS OF CUBIC POLY- 
NOMIAL 


NUMBER OF DIFFERENT REAL ROOTS 
ARRAY CONTAINING REAL ROOTS 


01834: 

IF<TE8T .LT »0) THEN 

01833: 

NROOTS-3 

01836: 

PHI-ACOS < -SIGNBtSQRT ( BB/ < -AAA > ) ) 

01837: 

SRT»SQRT(-A/3.) 

01838: 

RR ( 1 ) -2 . tSRTfCOS < PHI /3 . ) 

01839: 

RR<2)-2.*SRT*C08<PHI/3.+2.tPI/3.) 

01840: 

RR ( 3 ) *2 ♦ HSRTtCOS < PHI /3 ♦ +4 . *PI/3 # ) 

01841: 

ELSE IF<TEST .QT .0) THEN 

01842: 

NROOTS-1 

01843: 

SI — ♦ 5IB+SQRT ( TEST > 

01844: 

S2-- • 5IB-SQRT ( TEST ) 

01843: 

IF <ABS(Sl).LT.1.0E-06> THEN 

01846: 

SIGNS1*0. 

01847: 

ELSE 

01848: 

SIGNS1-P1/. ;BS<51 ) 

01849: 

ENDIF 

oiaso: 

IF (ABS<S2).LT.1.0E-06) THEN 

01831: 

SIGN32-0. 

01852: 

ELSE 

01853: 

SIGNS2*S2/CIBS<S2) 

01834: 

ENDIF 

01853: 

RR<l)=SIGNSlt<ABS(Sl>**< 1./3.) ) 

01856: 

* +SlGNS2*(ABS(S2)**<l./3. >) 

01857: 

ELSE 

01858: 

NR00TS-2 

01859: 

RR<1)>- SIGNB*2.*SQRT<-A/3. ) 

01860: 

RR(2)=SIGNB*S0RT(-A/3.) 

01861: 

END IF 

01862: 

DO 10 1-1,3 

01863: 

RR( I )=RR< I )-P/3. 

01864: 

10 CONTINUE 

01865: 

RETURN 

01866: 

01867:c 

01868:c 

0186 ?:c 

END 

01870: 

0187i:c 

01872:c 

SUBROUTINE EXTREMES<X r Y,TMiiX,TMIN,NR,NC> 

01873:c 

EXTREMES FINDS THE MAXIMUM AND MINIMUM VALUES 

01874JC 

01875:c 

AMONG THE ELEMENTS OF TUO THO-DIMENSIONAL ARRAYS 

01876:c 

018771C 

INPUT 

018785C 

X...X COMPONENT 

01879JC 

Y...Y COMPONENT 

01880 :c 

NR... DIMENSION (JF X ARRAY 

0188KC 

01882 :c 

NC... DIMENSION OF Y ARRAY 

01883 : c 
01884:c 

OUTPUT 


01885.C 


TMAX... MAXIMUM VALUE IN ARRAYS 

01886.’C 


THIN... MINIMUM VALUE IN ARRAYS 

018871C 



01888: 


REAL X( 100; 100) fYdOOf 100) 

01889: 


TMAX*X(1;1) 

01890: 


THIN*X(1,1) 

01891: 


DO 20 I*1»NR 

01892: 


DO 10 J* 1 t NO 

01893: 


IF ( X ( I f J ) . GT . TMAX ) TMAX*X<I,J) 

01894: 


IF<X(I»J)»LT .THIN) TMIN*X<I,J) 

01895: 


IF ( Y ( I . J ) . GT . TMAX ) TMAX*Y<I,J) 

01896: 


IF(Y( I f J) .LT.TMIN) THIN*Y(I,J) 

01897: 

10 

CONTINUE 

01898: 

20 

CONTINUE 

01899: 


RETURN 

01900: 


END 

oi9oi :c 



01902:c 



01903:c 



01904: 


SUBROUTINE NORM ( X , Y , TMAX , TMI N . NR , NO ) 

0t905:c 



01906:c 

THIS ROUTINE NORMALIZES THE VALUES OF TWO 

01907:c 

TWO 

DIMENSIONAL ARRAYS SO THAT THEY LIE 

01908:c 

BETWEEN 0 AND 1 INCLUSIVE. 

01909:c 



oi9io:c 

INPUT 

01911JC 



019121C 


X...X COMPONENT ARRAY ON INPUT AND 

01913JC 


NORMALIZED X COMPONENT ARRAY ON OUTPUT 

01914JC 


Y...Y COMPONENT ARRAY ON INPUT AND 

01915:0 


NORMALIZED Y COMPONENT ARRAY ON OUTPUT 

019161C 


NR... DIMENSION OF X ARRAY 

01917:c 


NC... DIMENSION OF Y ARRAY 

0191810 


TMAX... MAXIMUM VALUE IN ARRAYS 

01919:0 


TMIN... MINIMUM VALUE IN ARRAYS 

01920*0 



01921 1 C 

OUTPUT 

01922JC 



0192310 


X... NORMAL I ZED X ARRAY 

01924 :c 


Y... NORMAL I ZED Y ARRAY 

01925:c 



01926*. 


REAL X( 100. 100) »Y( 100; 100> 

01927J 


DO 20 1*1. NR 

01928: 


DO 10 J=1,NC 

01929: 


X(I,J)«(X(I.J)-TMIN)/<TMAX-TMIN) 

01930: 


Y< I. J)*< Y ( I , J)-TMIN)/(TMAX-TMIN) 

01931: 

10 

CONTINUE 

01932: 

20 

CONTINUE 

01933: 


RETURN 

01934: 


END 

BOTTOM 
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